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On Trails of Axisymmetric Hypersonic Blunt 
Bodies Flying Through the Atmosphere 


SAUL FELDMAN* 
Electro-Optical Systems, Inc. 


Summary 


The trail left in the atmosphere by a body moving at hypersonic 
speeds is the subject of theoretical treatment. The times re- 
quired for ionization and dissociation (and their inverse processes) 
to go to completion, when compared to the flow times of a gas 
particle, are important in determining the observable effects of 
hypersonic trails—i.e., emitted thermal radiation and reflection 
of electromagnetic waves from the trail. 

In order to simplify the theoretical treatment, the trail is 
divided into two regions: (1) the expansion-controlled trail, 
which treats the behavior of the wake behind the body up to a 
point, along the direction of flight, where the pressure decays to 
the free-stream value and cooling is controlled principally by the 
expansion of the flow, and (2) the conduction-controlled trail, 
where the trail cools mainly by diffusion of he#t away from the 
high-temperature core. 

The influence of the details of the body shape on the ob- 
servables are discussed and a simple computational procedure for 
the behavior of the conduction-controlled trail is developed based 
on integral methods. Results of calculations that assume thermo- 
dynamic equilibrium of the flow field give the values of the 
thermodynamic variables in the trail of a sphere, axial distribu- 
tions of emitted thermal radiation, and maps of electron density 
distribution. It is shown that the cooling of the conduction- 
controlled trail is essentially due to conduction of heat and that 
viscous effects are not important. It is found that this portion 
of the trail does not widen as one proceeds downstream. Flight 
velocities considered vary between 15,000 and 35,000 ft/sec and 
altitudes range between 100,000 and 250,000 ft. 


Symbols 
d += quantity defined in Fig. 2 
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= u/u.., defined after Eq. (6) 
= heat flux 
total radiation intensity per unit length of trail (Eq. 4) 
three-body recombination rate constant 
length of frozen trail 
= recombination length 
pressure 
radius 
r/¥n 
temperature 
time 
flight velocity 
= length along the trail measured from the center of the 
spherical nose 
= x,/Tn 
shock detachment distance 
boundary-layer thickness 
quantity defined in Fig. 2 
emissivity of air per unit thickness 
density 
= angles defined in Fig. 9 


Subscripts 
= free-stream conditions 
density conditions at 0°C and 1 atmosphere pressure 
= denotes conditions in the expansion-controlled trail 
denotes conditions in the conduction-controlled trail 


= nose 
(1) Introduction 


6 ew INTEREST IN hypersonic flow about blunt bodies 
has been stimulated in the last few years by 
intercontinental ballistic missile and satellite programs. 
Meteor physicists have been interested for many years 
in the trails left by meteors.’ ? The meteor inter- 
action with the atmosphere occurs at sufficiently high 
altitudes so that free molecular flow is of importance. 
Recently, interest in the trails left in the atmosphere by 
a body moving at hypersonic speed has also been 
stimulated by the desire to monitor man-made objects 
entering from space. In this second case, owing to the 
lower flight velocity than in the case of meteors and to 
the larger size of the body, the important interaction 
occurs in the continuum flow regime. 
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Among the effects that can be observed during 
atmospheric entry, either in the case of meteors, 
meteorites, or man-made bodies, are the thermal radia- 
tion emitted by the hot gas left in the trail and the re- 
flection by the ionized material of man-generated 
electromagnetic radiation. For the purpose of making 


the problem amenable to a simple analysis, we will 
neglect here the effect that ablation and injection of 
materials into the flow field has on the observables. 

When the flight regime is such that a shock exists in 
front of the body, high thermal and kinetic energy is 
imparted only to the gas that flows through the blunt 
portion of the shock. At high altitudes most of the gas 
that goes through the blunt portion of the shock, re- 
gardless of body geometry, flows into the body bound- 
ary layer. Under these conditions, a good portion of 
the energy will be transferred to the body and the rest 
will have to be dissipated in the trail. At lower al- 
titudes, the boundary-layer thickness becomes much 
smaller than, for example, the shock detachment dis- 
tance. Under those conditions, most of the energy 
is contained in the inviscid portion of the flow 
field, and the effects of the body viscous boundary 
layer on the features of the trail could be neglected. 
However, if the blunt portion of the body has a long 
slender afterbody, the boundary layer on the slender 
portion may entrain, even at low altitudes, most of the 
high-energy gas, thus leaving little energy to be dis- 
sipated in the trail due to pressure drag; however, 
friction drag may then become of importance. The 
problem to be investigated here will be restricted to rel- 
atively short blunt bodies,—i.e., where the boundary- 
layer growth on the afterbody has a small effect on the 
amount of energy left to be dissipated in the trail. 

In order to design adequate detectors of either the 
thermal radiation emitted by the gas in the trail, or its 
electromagnetic reflection properties, an estimate is 
required of the temperature and density field in the 
trail. Before *attacking this problem in detail, it will 
be of interest to point out the physicochemical phe- 
nomena that may arise. The observables will strongly 
depend on the chemical composition of the gas left 
behind in trail. The gas chemistry will be influenced 
by the relative magnitude of the gas particle flow times 
and the times necessary for processes like molecular 
vibration, dissociation, ionization, and their inverse 
processes to go to completion. Once the chemical 
composition of the gas flowing off the body and into the 
trail is established, the problem is to determine what 
portion of the trail will be laminar and what portion 
turbulent. This is a very important question, since it 
is well known that transfer processes in turbulent flow 
at high Reynolds numbers are thousands of times more 
rapid than in laminar flow. Hence, a high-tempera- 
ture turbulent wake will probably cool faster than a 
laminar one by several orders of magnitude, resulting 
in a shorter observable trail than if it were laminar. 

Fig. 1 is the result of an attempt to summarize some 
estimates of the physicochemical processes of impor- 
tance in determining the observables in a trail in the 
continuum gas flow regime. We shall first point out 
the meaning of the lines and later discuss their implica- 
tion. Using data obtained from chemical kinetic ex- 
periments in air carried out at our Laboratory, Teare* 


* J. D. Teare of our laboratory has kindly made available his 


results to the author. 
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has calculated the distance d (Fig. 2) behind a normal 
shock necessary for the gas temperature to relax down 
to a value that is 35 percent above the equilibrium 
value. In order for the temperature to drop to within 
10 percent of the equilibrium value, a distance of the 
order of 10d would be needed. Lines of constant value 
of d have been plotted in Fig. 1; at sufficiently high 
altitudes the dissociation distance d becomes larger 
than the shock detachment distance 6 which is, to first 
approximation, equal to r,/10 where 7, denotes the body 
nose radius. At sufficiently high altitudes, d/6 might 
be, for situations of interest, larger than unity; thus 
the flow over the nose would be out of equilibrium. 
However, for large meteors the flow will be in equilibrium 
at those altitudes. Similarly, lines of constant e, 
the distance needed to reach 120 percent! of equi- 
librium electron density », (from below, cf. Fig. 2), 
have also been included in Fig. 1. 

In the attainment of equilibrium in the forward 
portion of a body, two-body collisions are involved. 
Since recombination in the expansion region of the 
flow is a three-body process, it is possible for the flow, 
although in equilibrium in the compression regions of 
the flow, to be out of equilibrium in the expansion regions 
of the flow field. Thus, constant values of the trail length 
/,, where the flow is frozen, have been estimated (cf. 
Appendix A) and plotted in Fig. 1. 

Once the gas flow expands around the body into the 
trail and cools, it begins to recombine. The energy 
released during atomic recombination tends to in- 
crease the temperature of the trail’s core, and heat 
conduction radially away tends to cool the core. By 
use of a recombination rate constant k, (Appendix) 
of 10% [4,500/7 (°K)'/? mole? cm~* sec~!, an estimate 
can be made of the recombination heating in the trail, 
as well as the rate of radial heat conduction. These 
rates will depend respectively on the cube of the gas 
density (three-body collisions) and on the nose radius, 
because of the radial temperature gradient. Equat- 
ing both rates leads to an expression for the free- 
stream density p./po(pp = 1.29 X 107% g/cm’) 
i.€., approximately 


Po/ Py = 2 ~*~ 10 2 [r,(cm) | -1/3 (1) 


Thus, the altitude at which the heat of recombination is 
balanced by heat conduction is fairly insensitive to nose 
radius, since it depends on r,'/*. The dependence of this 
altitude on recombination rate constant is also through 
a cubic root. In Fig. 1 a line has been drawn above 
which, for a 10-cm nose radius, the heat fluxes q sat- 
isfy the relationship Grecomb. O.1dconduct. At altitudes 
above that line, if the trail is in the laminar flow regime, 
it will cool to ambient temperature before recombining. 
Trajectories of typical IRBM and ICBM nose cones 
as well as a satellite trajectory are alsoshown. We will 
now discuss the implications of Fig. 1. 

For bodies of usual size and at altitudes above 250, 


+ The fact that different overshoots were chosen for d and e 
is immaterial for the present discussion. 
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frozen flow. 


000 ft. d and the recombination length /, may be larger 
than the physical size of the body of interest: thus the 
flow over large parts of the body will be out of equi- 
librium. This conclusion holds almost independently 
of flight velocity. At lower altitudes the flow over the 
nose will be in equilibrium; however, the flow in the 
trail will be chemically frozen over distances of the 
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order of thousands ef feet. In order to get an under- 
standing of the variation of properties of the gas in the 
trail, which is between the extremes of chemical equi- 
librium and frozen states, we will next refer to Figs. 
3-7. Figs. 3-5 give the values of the thermodynamic 
coordinates, if they were in thermodynamic equilibrium, 
of the gas on a streamline, that has gone from free 
stream through a normal shock at the flight velocity 
and expanded isentropically around the body to free- 
This, of course, neglects variations 
To a first 


stream pressure. 
of entropy downstream of the nose shock. 
approximation, it can be seen from Fig. 3 that the 
temperature increases linearly with flight velocity. 
The temperatures along the axis of symmetry of the 
flow field, upstream of the point where the pressure 
decays to the free-stream value, will be higher than 
those indicated in Fig. 3, the highest temperature oc- 
curring at the forward stagnation point of the blunt 
body. The stagnation point temperatures’ are between 
2,500 and 5,000°K higher, depending on velocity, than 
the ones shown in Fig. 3. Fig. 4 gives the ratio of the 
density to the free-stream density, and it can be seen 
that at low velocity, almost independent of altitude, 
the density is about 1/10 the free-stream value, while 
at higher velocities it goes down to about 1/30 of free- 
stream density. The velocity of the flow (Fig. 5), 
with respect to a reference frame fixed in the body, is 
about 75 percent of the flight velocity regardless of its 
absolute value. Figs. 6 and 7 compare, at an altitude 
of 200,000 ft, the temperatures and densities that are 
obtained by assuming equilibrium flow across the nor- 
mal shock and frozen or equilibrium flow in the expan- 
sion; altitude variation has a small effect on the relative 
differences shown. In the frozen case the gas was as- 
sumed to freeze (vibrationally and chemically) at two 
different points on a sphere; in one case it was frozen 
near the sonic point while in the other it was frozen 
at 90° from the stagnation point. The reason for 
choosing two different values is that the gas will 
freeze chemically at different locations depending on the 
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body size. If the flow across the shock is in equilibrium 
but later freezes when it expands around the body into 
the trail, the temperatures can be radically different 
from the equilibrium value (Fig. 6); this will have a 
strong effect on the trail observables. From Fig. 6 
it can be concluded that if the flow is allowed to re- 
combine at free-stream pressure, the temperatures could 
increase by a factor which will vary between 2 and 5 
depending on the flight velocity. Figs. 6 and 7 also 
have lines designated as y = 1.4. They denote the 
state of the gas if it were vibrationally and chemically 
frozen across the nose shock and throughout the flow 
field; in that case it is interesting to note that the tem- 
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perature and density come quite close to the equilibrium 

value. However, the gas composition will be radically 
different and, if left to itself, the y = 1.4 gas will dis- 
sociate at constant pressure as it flows in the trail, lead- 
ing to a decrease in temperature. 

Returning to Fig. 1, and with the background just 
presented, we can see that at high velocities and al- 
titudes the trail will probably be cold (frozen) and lam- 
inar, and it will probably become turbulent before un- 
freezing. Since above altitudes of the order of 110,000 
ft heat is lost faster from the frozen wake than it is 
generated by atomic recombination, the trail will 
reach free-stream temperature before the recombina- 
tion process will rewarm it. At sufficiently low al- 
titudes, the flow will be in thermodynamic equilibrium 
and the trail will first cool due to expansion from the high 
pressure in the neighborhood of the body to the free- 
stream value, at that point reaching the temperatures 


given in Fig. 3. From that point on (denoted as the 
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pb > Pp. point), the trail will cool by laminar or tur- 
bulent convection. 

Summarizing, if only laminar flow is considered in 
regard to the equilibrium trail, its temperature will 
monotonically decrease with distance downstream from 
the body. In the case of the frozen wake, the tempera- 
ture will first be lower than the equilibrium tempera- 
ture and further downstream, depending on the body 
size and altitude, the wake could either reheat and then 
cool again or the frozen chemical species could just 
diffuse away without appreciable reheating. If in these 
considerations we include the fact that the trail some- 
where along its path could become turbulent (thus 
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Fic. 10. Radiation intensity of gas in region 3 as a function of 
flare angle (4), and expansion angle (6; + 42). 


increasing drastically the convective process), the com- 
plete elimination of the reheat process of the frozen 
wake could result. The local chemical composition of 
the trail, together with the temperature and density, 
will determine the thermal radiation emitted by the gas. 

From the lines of constant values of d and e (Fig. 1), 
when the altitude is sufficiently large, the possibility 
exists of not reaching equilibrium at the stagnation 
point; thus the electron density in the entire flow field 
may be out of equilibrium. If the altitude is de- 
creased, equilibrium ionization is reached at the stagna- 
tion point, but from there on downstream the axial 
electron density distribution in the trail will depend on 
the magnitudes (if the flow is laminar) of the rates of 
ionic and atomic recombination times relative to the 
flow times. At high altitudes both recombination rates 
are low and thus the electron density in the trail stays 
high above the equilibrium value. At intermediate 
altitudes the ionic recombination (NOt + e > N + 
O) is greater than the pertinent atomic recombination. 


At these altitudes the electron density will be in equi- 
librium with the local frozen temperature, thus over- 
shooting the final value which would be obtained after 
the gas recombines (Fig. 8),* i.e., equilibrium. How- 
ever, if the flow were turbulent, the electron density 
might stay low and never recover. 

Conversely, if the ionic recombination is slow com- 
pared to atomic recombination—i.e., low altitudes and 
velocities—the electron density will always be higher 
than the equilibrium value (Fig. 8). Again the final 
fate of the trail is determined by ambipolar,* laminar, 
or turbulent diffusion. 

As can be seen from the foregoing discussion, the 
possible variations of the physicochemical phenomena 
in the trail are many. 

In order to develop a computational method for 
determining the flow field in the trail, to acquire an 
understanding of the scaling laws that govern the 
trail observables, and to estimate their order of mag- 
nitude, it will be useful to choose for investigation a 
particular case out of the multitude of possible situa- 
tions. We have chosen to investigate first a compara- 
tively simple case: the behavior of the laminar trail 
of axisymmetric, nonablating blunt bodies when the 
flow is in thermodynamic equilibrium. We will first 
treat some generalities of the inviscid flow field in the 
neighborhood of the body because it constitutes the 
initial condition for the treatment of the trail. We 


* The arguments leading to Fig. 8 were first given by S. C. 
Lin of our Laboratory in an unpublished paper. 
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laminar trail about a sphere. 


will next discuss the behavior of the laminar trail; 
thermal radiation and electron density distributions 
throughout the flow field of a typical blunt body will 
be presented.* 


(2) The Flow Field in the Neighborhood of a 
Blunt Body 


In the case of blunt-nosed bodies at sufficiently high 
speeds, most of the entropy increase of the gas occurs 
when a streamline goes across the normal part of the 
shock near the stagnation region. The specific details 
of the flow over the afterbody should cause only minor 
changes in the entropy level. 

(2.1) Effect of the Details of the Blunt-Body Geometry on 
the Observables 

The general appearance of man-made re-entry ob- 
jects will be as depicted in Fig. 9 where a general two- 
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dimensional, blunt-nosed body is shown. It may be 
short as indicated by the cross-hatched section, or 
long as shown by the shaded section. In general the 
body may have a flare denoted by the angle 4, that 
will further shock the gas that comes through the nose 
wave. (The necessity of the flare is dictated by aero- 
dynamic stability considerations.) The gas expands 
beyond the base of the body flowing by the dead water 
region, and then shocks again so that the direction of 
flow downstream of the last shock is parallel to the 
flight direction. If the flow process between region | 
(Fig. 9) and region 3 were isentropic, the body geom- 
etry would have little effect on the state of the gas'in 
region 3, which determines the initial state of the gas in 
the trail proper. However, the process is not isentropic 
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trail of a sphere. 


and it is of interest to determine, for example, how the 
state of the gas that has gone through the normal 
portion of the nose shock and arrived at region 3 is 
affected by details of body geometry. For that purpose 
we have considered the two-dimensional body of Fig. 9 
in which the gas is shocked by a wedge of angle 4; 
expanded by an angle 6, + 4, and further shocked by a 
wedge of angle @. The angle 6, has been varied between 
0 and 30° which encompasses most of the possibilities 
of practical interest, and 6. was varied between () and 
20° which essentially covers the expansion angles that 
have been found experimentally in wake flows reported 
in the literature (see for example reference 5). The 


* The behavior of the technically important chemically react- 
ing trail will be presented in a future paper. 
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calculation was carried out for a velocity of 20,000 
ft/sec and an altitude of 100,000 ft. The air was as- 
sumed to be in thermodynamic equilibrium.* Cal- 
culations were made for the case of a very short body, 
for which the pressure upstream of the wedge 6; was 
5 percent of the stagnation pressure—i.e., pi/Ps: = 
0.05—and for a long body where ~;/ps, = 0.01. The 
two-dimensional case where 6; = 6. = 0 and the flow 
is not deflected corresponds to the axisymmetric case 
of a hemisphere-cylinder. In order to evaluate the 
effect of the temperature and density variations due to 
geometry on thermal radiation in the trail, the equi- 
librium radiation emitted in the optical (22,000 
10,000 A) and infrared (5.3 u) regions of the spectrum 
has been calculated using the tables of Kivel and 
Bailey.6 The results are presented in Fig. 10. Once 
a value ~;/ps: is chosen (and this is essentially de- 
termined for a given body by the axial distance of the 
flare from the stagnation point), the radiation and the 
electron density will not vary by more than a factor of 
3. This means that any conclusions as to magnitude 
of thermal radiation and electron density that are 
arrived at on the basis of calculation for a simple blunt- 
body geometry like a hemisphere-cylinder, that does not 
include the geometrical details of the afterbody, may 
be off by a factor of not more than 3. Although this 
factor may be of importance in engineering applica- 
tions, it certainly can be disregarded in an attempt first 
to understand the order of magnitude of the observ- 
ables of interest. 

One other class of re-entry body geometry can be 
conceived to be composed of sphere-capped cones. 
When the cone is of small included angle and long rel- 
ative to the nose radius, the column of gas that is 
heated near the stagnation region will have been in- 
jected into the body boundary layer by the time it 
leaves the body. The study of such an effect on the ob- 
servables is beyond the scope of this paper. 


(2.2) Estimate of the Importance of Viscous Effects on 
the Initial State of the Trail 

Thus far, we have been discussing the inviscid flow 
field only. We will next examine the influence of 
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viscous effects on the foregoing results. For the pur- 
pose of estimating these effects, let us assume that the 
body of interest is a sphere. By use of the theoretical 
work of Kemp, ef al.’ for the equilibrium air laminar 
boundary layer away from a stagnation point, it is pos- 
sible to correlate the calculations for the boundary- 
layer thickness 6,; at the shoulder (90° from the stagna- 
tion point) of a sphere by the following formula: 

bn | 2 X 10-* Vem F 

', Vp. ‘po)%n (cm) *) 


By use of this formula, one can estimate the volume 
occupied by this boundary layer by the time it is ex- 
panded from the region of high pressure on the body to 
the free-stream value. To a first-order approximation, 
it is possible to use a one-dimensional treatment and 
assume no mixing in the dead-water region behind the 
body. The velocity of the flow changes only by a 
small amount when going from the shoulder to free- 
stream pressure, and the density decreases by a factor 
of about 10. Since the cross-sectional ring area of the 
boundary-layer flow at the shoulder becomes a cir- 
cular cross section when it comes together on the axis 
of symmetry, the following result follows from con- 
servation of mass 


(6, Yn) p > po _ V 20(6 Wd eee (3) 


Thus, if the boundary layer has a small effect on the 
inviscid flow of the forward portion of the body (which 
is true for a blunt body as defined in the introduction), 
it will also occupy a negligible volume once it is ex- 
panded to free-stream pressure. Thus, the boundary 
layer can be neglected when considering its effects on 
the observables. This will only be true if we neglect, 
as we will here, the effect on the observables of abla- 
tion and injection of materials into the flow field. 


(2.3) The Inviscid Flow Field About a Typical Blunt Body 


We have shown above that in the case of blunt bodies, 
the details of the body geometry can be neglected pro- 
vided account is taken of the effects of the strong bow 


wave ahead of the body. Thus, the flow about a sphere 
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could be considered to be representative of a typical 
blunt body. However, the calculation of the in- 
viscid flow about a sphere is more laborious (because of 
the trailing shock) than, for example, the flow about 
a long hemisphere-cylinder. The details of such a 
flow field at a flight velocity of 17,500 ft/sec and an 
altitude of 60,000 ft in equilibrium air has been treated 
by the author. The results presented in Figs. 2, 3, 
4, and 8 of reference 8 were used in the radiation es- 
timates made in the next section. 

Owing to the fact that the gas density near the body 
is very low (Fig. 3 of reference 8), one could deform 
the cylindrical portion of the body without affecting 
the flow field very much. Thus, for example, if one 
were to eliminate the cylindrical portion of the body 
altogether, the profiles of temperature at any given 
axial station far from the nose could be corrected merely 
by displacing the values on the body toward the axis 
of the flow. The entire temperature profile could then 
be obtained by fairing the outer portion to the newly 
determined value on the axis. A more exact correction 
could be obtained by accounting for the fact that the 
flow area changes when the cylindrical portion of the 
body is eliminated. This effect will lead to a slightly 
lower temperature on the axis than the one obtained 
from the approximate method just mentioned. For 
our purpose, it will be acceptable to ignore this small 
temperature change. 


(2.4) Thermal Radiation and Electron Density Distribu- 
tion in the Expansion-Controlled Trail of a Blunt Body 

In order to determine the temperature and density 
in the flow field, calculations of the flow field (by the 
method described in Section 2.3 and in reference 8) 
have been carried out at altitudes of 100,000 and 250, 
000 ft and at flight velocities varying between 15,000 
and 25,000 ft/sec. 

With a knowledge of the temperature and density 
distribution in the trail, it is possible with the help 
of the data given by Kivel and Bailey® to determine the 
emissivity per unit thickness, «/Z, of high-temperature 
air. Since under the flight conditions of interest air 
is transparent, the emissivity will turn out to be much 
less than unity. The thermal radiation intensity per 
unit length of the trail, Z;, can be determined from the 
following equation: 


ijt? = tno f : T‘*r* dr* (4) 
es 


where o is the Stephan-Boltzmann constant given as 
o = 5.672 X 10~'” watts/cm?-°K‘4, r* denotes the ratio 
of the radial coordinate 7 (normal to the axis) to the 
nose radius 7,, and the upper limit of integration has 
been changed from the shock-wave coordinate to in- 
finity. It should be noted that the radiation intensity 
per unit length, /;, is directly proportional to the cross- 
sectional area of the high-temperature gas column 
and therefore proportional to 7,” because of the relation- 
ship between 7, and the radius of the hot gas column; 
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Fic. 16. Temperature along the axis of an infinite cylinder of 
infinite radius cooled by heat conduction. 


Kivel and Bailey® give the emissivity per unit length 
integrated over the optical spectrum (2,000—10,000 
A) and also the emissivity per unit length in the in- 
frared region (5.3 uw) due to vibrational and rotational 
transitions of the nitric oxide molecule (NO). Fig. 
11 shows the radial distribution of the radiation in- 
tensity in the infrared range for a typical case. Itisa 
plot of the integrand of Eq. (4) for a sphere in the ex- 
pansion-controlled trail at a station 80 nose-radii away 
from the shoulder; the flight velocity is 20,000 ft/sec 
and the altitude 100,000 ft. The area under the curve 
gives [ino/?rn?. The integrand is small near r* = 0) 
because 7* appears as a factor in Eq. (4). This means 
that errors made in the calculation of the intensity, 
due to the existence of viscous effects near the axis 
that are remnant from the body boundary layer, will 
be decreased. The results of finding the area under 
many curves similar to the one given in Fig. 11 are 
present in Fig. 12 which gives the radiation intensity 
I;, per cm length of trail, normalized with respect to 
the nose radius 7, (cm) squared, as a function of the 
dimensionless length x,* = x;/r,; x,* is assumed to 
be zero at the shoulder of the sphere. 

Thus the absolute radiation per unit length will be 
proportional to the square of the nose radius. The 
total integrated radiation in the trail up to the point 
where the pressure has decayed to the free-stream 
value (we will denote this point along the axis as the 
pb — p,, point) will be proportional to the cube of the 
nose radius. Fig. 13 gives the total integrated radia- 
tion as a function of flight velocity. The total radia- 
tion over the whole spectrum from the gas in the nose 
cap up to the shoulder of the sphere has also been in- 
cluded in Fig. 13. It was obtained by multiplying the 
stagnation-point radiation per unit volume, given by 
Kivel and Bailey, by an effective volume of 10~'r,’. 
(This effective volume was found to give roughly the 
correct amount of total radiation from a nose cap.) 

Equilibrium electron densities can be obtained by 
knowing the temperature-density field; Figs. 14 and 
15 show the results. Ambipolar diffusion of elec- 
trons will be important only in determining the elec- 
tron maps for small bodies at high altitudes (200,000 
ft and above). This effect is neglected here. 
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(3) The Conduction-Controlled Laminar Trail 


Once the pressure along the flow direction has de- 
cayed to the free-stream value, further cooling of the 
wake will be governed by heat conduction from the 
(Radia- 
tion cooling can be shown to be negligible.) We have 
thus far neglected heat transfer in the expansion region. 
The justification for this is that at altitudes not much 
below 200,000 ft in laminar flow, heat transfer will 


high-temperature core to the outer cold fluid. 


usually be small compared to the expansion cooling, 
and thus will be unimportant in determining the prop- 
erties of the expansion-controlled trail. The observa- 
ble effects in general will be small at high altitudes. 
As an example (cf. Section 3.4), the heat transfer and 
expansion cooling will have to be calculated simul- 
taneously* at an altitude of 250,000 ft and a velocity 
of 25,000 ft/sec when the body is smaller than 3 or 30 
em, depending on the observable of interest. The 
flight conditions for which cooling due to heat conduc- 
tion and expansion are of the same order of magnitude 
will be determined automatically from the calculations 
of the conduction-controlled trail—i.e., when the con- 
duction-controlled trail cools in distances of the same 
order as the length of the expansion-controlled trail. 
In this paper we will neglect this interaction since a 
correction can always be applied on the basis of the 
results of our calculations. For the purpose of develop- 
ing a method for calculating the cooling of the con- 
duction-controlled trail, we will first limit ourselves to 
a study of the behavior of a trail whose initial condi- 
tion—1.e., the properties at the point p > p«—are de- 
termined from the characteristic calculation referred 
to before; the hemisphere-cylinder flying at 17,500 
ft/see and 60,000 ft altitude. 


(3.1) The Cooling of an Infinite Cylinder with Gaus- 
sian Temperature Distribution and Constant Physical 
Properties 

Since we have a case of heat convection in a 
compressible fluid with variable properties, we would 


* The necessity of treating the combined effect at high altitudes 
was first pointed out by R. J. Goulard. 
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first like to estimate the effect of the variable prop- 
erties on the cooling rates. For this purpose let 
us ignore the fact that we are dealing with a shear 
flow. Let us assume two different 
thermal diffusivity of the material: one determined 


values for the 


by the properties on the axis of the characteristic cal- 
culation, and the other determined by the conditions 
far from the axis—i.e., free-stream conditions—this 
should determine the range within which the cor- 
rect solution would fall. Naturally, owing to the high 
temperature on the axis, the diffusivity will be higher 
there and thus cooling will be more rapid than when the 
Time ¢ in the diffusion 
In this case the 


free-stream diffusivity is used. 
problem can be replaced by x/u 
results can be correlated as a function of x/r,? instead 
of the usual //r,*. Assuming a Gaussian radial dis- 
tribution of temperature and the Gaussian width o 
to be one nose radius, the temperature variation along 
the axis of the cylinder can be found.’ The result is 
shown by (a) and (b) of Fig. 16. It can be seen that 
there is quite a range of uncertainty left by this simple- 
minded approach. In order to decrease the un- 
certainty, it will be necessary to include the variable 


properties of the flow field. 


(3.2) The Cooling of an Infinite Cylinder with Gaussian 
Temperature Distribution and Variable Properties 

There are two ways in which the variable properties 
in the cooling of a cylinder can be accounted for. One 
of them involves a full numerical integration of the 
heat conduction equation; this is quite a laborious 
The other approach is to use integral tech- 
This last approach has been used here and the 


method. 
niques. 

results for the present case, for a Prandtl number Pr 
of 0.71 and a Gaussian width in a Howarth transformed 
coordinate system [Eq. (8)] of one nose radius, is 
given by curve (c) in Fig. 16. 

In order to investigate the importance of the dis- 
crepancy between the results given by curve (c) of 
Fig. 16 and the correct solution, the full problem in- 
cluding the convection and viscous terms will have to 
be investigated (this was done in reference 10, and 
will be shown in Section 3.3). However, since, as 
it will be seen, the errors involved in curve (c), Fig. 
16, are relatively small—i.e., less than 5 percent in tem- 
perature at a given axial station—it will be of interest 
to summarize here the approach used in obtaining 
curve (c). It should be mentioned at this point that 
according to Kivel and Bailey,® a 10 percent error in 
temperature leads to a 100 percent error in the radiant 
intensity. 

The cooling of an axisymmetric column of gas with 
variable properties is governed by the following differ- 
ential equation: 

10 por or . 


= Co — (90) 
r or’ or me Ot 


where r denotes the radial coordinate, h the thermal 
conductivity, 7 the temperature, p the density, cp 
the heat capacity at constant pressure, and ¢ the time. 
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The 0/o¢ in Eq. (5) could in our case be replaced by 
u(0/Ox) where a is some average velocity of the hot 
gas column with respect to the moving vehicle and x is 
the axial coordinate. For hypersonic flight this ve- 
locity could be taken (Fig. 5) to be about 75 percent of 
the free-stream velocity. Since the characteristic 
dimension of the heated gas column is determined by 
the nose radius 7,, it is convenient in Eq. (5) to let 
r* = r/r,. Rewriting (5) in terms of static enthalpy 
h and Prandtl number Pr = ywc,/k, where » denotes 
the viscosity coefficient (after referring some quantities 
to free-stream values) leads to 


1 0 E m ee p . O(h/h.) (6) 

La ‘Lis | [El Ad Adda tM 

v* Or* LPr uo or* p. O(x*/Re) 
where x* = x/r,, Re is a Reynolds number based on 
free-stream conditions (Ke = pafa%,/Mo), and f is 


a factor which is in the neighborhood of 0.75 depend- 
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Fic. 18. Gaussian fit of the enthalpy profile at the p > po 
point of the characteristic solution. 


ing weakly on flight velocity and altitude (Fig. 5). 
The function of the Reynolds number here is not to 
account for viscous effects but to attempt to extract 
the density dependence of the cooling distance [ef. 
Eq. (8) below]. 

Eq. (6) is identical to the relation that is obtained if, 
in the energy equation usually used in fluid dynamics, 
the radial velocity together with the radial gradient 
of the axial velocity are set equal to zero. It should be 
noted that, for a given flight velocity, p/p.. appearing 
on the right side of Eq. (6) has a weak dependence on 
altitude (Fig. 4). Eq. (6) will now be solved by use of 
an integral method; x will be assumed to be equal to 
zero at the initial value of the axial coordinate, which 
is taken as the point where the pressure on the axis, 
given by the characteristic calculation, first reaches 
the free-stream value p.,. (or any other prescribed value 
nearly equal to p..). This point will be denoted as the 
pb — p. point. Integrating Eq. (6) with respect to 


* 


r* from r* = 0 to r* = © and assuming that when r 
= 0 and r* = ©, all partials with respect to 7* are 


f y« ? O(h/h..) - (7) 
0 p. O(x*/Re) 


Using the Howarth transformation 


zero, leads to 


; \ (8) 
and letting | 
X = x*/Re 
and inserting them into Eq. (7) leads to 
a) 
; RdR=0 (9) 
OX Jo h 
From the density profiles at 1,* = 56 given in Fig. 3 


of reference 8, it is possible to obtain for a particular 
case the dependence of R upon r* by numerical inte- 
gration. The result is given in Fig. 17 as the curve 
x/r,” = 0. Let us assume next that the enthalpy is 
a Gaussian function of R given by 


(h/ho) — 1 = a(X) exp [—R?/B(X)] (10) 


Again from the profile at x* = 56 of Fig. 3 of reference 
8, it is possible to plot the enthalpy versus Rk. Fig. 
18 is a comparison between the profiles obtained from 
the characteristic calculation and several Gaussian pro- 
files of different widths. Since the important cooling 
will be determined by the enthalpy gradients close to the 
core, a Gaussian enthalpy profile could always be found 
that will fit the initial condition quite well. In addi- 
tion to the requirement of adequately fitting the in- 
itial data, the reason for the choice of a Gaussian 
representation for the enthalpy is that for large values 
of X the problem will reduce asymptotically to such a 
distribution. Also, and this is not the least reason, the 
Gaussian leads to an analytically integrable function.t 
Insertion of Eq. (10) into (9) yields 


a(X)B(X)/a(0)B(0) = 1 (11) 


where a(Q) and B(0) are obtained from the fitted re- 
sults just discussed. In order to solve for a(X) and 
@(X) in Eq. (11), one other relationship between them 
is required. This can be obtained by evaluating Eqs. 
(10) and (6) on the axis of symmetry—i.e., taking the 
limit of Eq. (6) asr* > 0. The result, after using Eq 
(11), is 


da 4 ua (X, 0) 


— = — —— ; 2(X) (12) 
dX Prfa(0)B(0) be " 


A rigorous solution of Eq. (12) for a(X) can be ob- 
tained numerically only. It will first be necessary to 
find the enthalpy, at each step, from Eq. (10). Know- 
t It should be noted, however, that because of the inadequate 
fit away from the core at the initial station, there is an inherent 
limitation as to how low an enthalpy can be reached in the 
present calculation and meaningful results still be obtained. 
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ing the enthalpy and the free-stream pressure, it is pos- 
sible to find the temperature and density of the gas at 
each point of the flow field by using the equation of state 
for air ora Mollier diagram.* The viscosity coefficient 
u can be determined from Sutherland’s formula as 


T \32 386 
© ae ( ) aes (13) 
Mo zis 113 + 7 


where 7’is the temperature in °K and yo is the viscosity 
of air at 273°K—.e, wo = 1.72 X 10-4 g/cm-sec. 

An approximate solution of Eq. (12) will be obtained 
next. The Prandtl number of air in thermodynamic 
equilibrium is approximately constant and equal to ().7. 
For the interesting range of variation of X, u/u will 
be a slowly varying function of X and can be assumed 
constant; we will denote it by the symbol Z@/y. Inte- 
gration of Eq. (12) under these circumstances leads to 


(X) 4 a(X, 0) t 
SS 2 E he sect Oem x (14) 
a(0) PrfB(0) 7 


Thus, Eqs. (10), (11), (13), and (14), together with 
the equation of state for air, solve the cooling problem 
completely. 


(3.3) The Entropy Wake Including the Convective and 
Viscous Terms 


In order to obtain a more adequate account of the 
observable phenomena in the trail of a hypersonic blunt 
body, it will be necessary to include in the calculation 
the convective and viscous terms that occur in the 
conservation equations. However, we are still not in- 
cluding the effect that the body boundary layer has 
when it is swept off the body and into the core of the 
trail. Thus the present calculations will, for example, 
at each axial station lead to an upper bound for the 
temperature. 

The interested reader can find a detailed description 
of the solution of the coupled momentum and energy 
equations in reference 10. The case of general value 
of the Prandtl number, as well as unity Prandtl number, 
is treated there. We will recall here some of the 
results obtained therein. The state and conservation 
equations were reduced to the following forms [Eqs. 
(15)—(19) |: 


<p = pZRT (state) (15) 

/ 7(p,h) and Z = 2Z(p,h) (Mollier chart) (16) 
b(X) = D/[a(xX)] {1 — [a(X)/2]} 

(momentum integral equation) (17) 


where a and 6 were the parameters in the assumed 
velocity distribution 


1 — (u/u.) = a(X) exp [—R?/b(X)] 


and D is a constant of integration essentially represent- 


ing half the drag coefficient of the body based on the 
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maximum cross-sectional area of the sphere constructed 
with the nose radius; 


d (a8 i as 4u..” (*\ in p(X, R)u(X, R) 
aU Uh he re 


R 
2 \ p >. 
} exp [—2R?/b(X yy} LI, o(X, =| dk? 


(energy integral equation) (18) 


and 
da 4 pw(X, 0) al(X) l 
dX Pr wpe B(X)1—a(X) 


(energy differential equation on the axis) (19) 


It should be noted that for the case where the viscous 
term is unimportant the right-hand side of Eq. (18) 
vanishes and Eqs. (18) and (19) then reduce, for the 
case of constant axial velocity, to Eqs. (11) and (14) 
of Section 3.2. Note that f of that section is replaced 
here by (1 — a), which is a constant if the velocity is 
constant. 

The right-hand side of Eq. (18), when evaluated for 
the flight case used here as a reference (alt. 60,000 ft; 
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u.. = 17,500 ft/sec), turns out to be of the order of 5 
percent of the other terms in the equation. As a 
numerical example, the conservation equations were 
solved in an electronic computer using Pr = 0.71. The 
temperature along the axis of the trail is given in Fig. 
19 by curve (a); the notation aB/a(0)B(0) # 1 with 
curve (a) indicates that this is the case when one solves 
the complete set of conservation equations. The dashed 
curve of Fig. 19 indicates the results obtained when the 
viscous effects are negligible—i.e., with a8/a(0)8(0) = 1 
replacing Eq. (18) and f assumed to be 0.8. Curve (c) 
of Fig. 16 is also included in Fig. 19 for comparison, 
from which it can be seen that it gives quite a good ap- 
proximation to curve (a). The inclusion of the viscous 
dissipation terms in curve (a) should and does give a 
lower temperature, at a given axial station, than if 
Varying the Prandtl number should 


they were left out. 
As an example, the 


have a small effect on the results. 
momentum and energy equations were decoupled by 
assuming Pr = 1 and then the momentum equation 
was solved by itself with a Gaussian width for the ve- 
locity profile of 0.75. The result coincides with the 
curve (c) in Fig. 19. 

Thus far we have presented only the temperature 
along the axis of the trail. Since we have finally ob- 
tained a fairly simple method of calculating the flow 
field, we will now present in Figs. 20 through 24 the 
profiles of the thermodynamic variables and_ their 
variation along the trail. These have been computed by 
use of Eq. (18). The dashed curves in Figs. 23 and 24 
are the results obtained when using a8/a(0)8(0) = 1. 

The reason for using (x2/7,”) (em~!) as a parameter 
is that it can readily be converted into the number of 
nose radii by multiplying the parameter by the nose 
radius incm. In Section 3.4 we will present some of the 
results for the observable quantities in terms of the 
general parameter X and also in terms of actual physical 
length. 

There is one striking feature that characterizes Figs. 
20 through 24, and it is that the profiles do not spread 
in the radial direction, as they usually do in diffusion 
problems. The reason can be found in the fact that 
the density is low in the core of the trail and is high away 
from it. The high-density outer region of the trail is a 
very efficient heat sink; in other words, a relatively 
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large decrease in temperature of the low-density core 
can be balanced by a small increase in temperature of 
the surrounding gas. Thus, the enthalpy is kept from 
spreading. Of course, this feature of the flow is some 
what exaggerated by the present calculation because, 
in the original fit of the enthalpy profile at x.,/7,2 = 0, 
as r* increased the Gaussian profile approached a 
smaller value of enthalpy than the true distribution 
(Fig. 18). Owing to the inadequacy of the fit to the 
original profile in the cool outer region, the present cal- 
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Fic. 22. 


culation becomes inaccurate when the temperature in 
the core decreases by more than 1,500°K. The density 
ratio p/p) (Fig. 24) does not approach asymptotically 
the free-stream value for the same reason. The fact 
that the Prandtl number used in the calculation was 
less than unity implies that energy diffuses faster than 
vorticity. This means that although at x. = 0 and « 
the total enthalpy is constant and equal to the free- 
stream value, for intermediate values of x2 it will first 
have to diminish in the core and then increase in the 
outer portions of the trail, and eventually for suffi- 
ciently large values of x2 the trend will have to reverse. 
Fig. 25 shows the total enthalpy profiles on an enlarged 
scale, in which the magnitude of the error at x. = 0, 
due to imperfect matching of the Gaussian profiles, can 
be seen. It should be noted that the profiles behave as 
expected and that for x2/r,” = 25,000 cm~', the trend 
on the axis of symmetry has reversed and the enthalpy 
there is already starting to increase. 

The percentage error introduced by using the integral 
method instead of integrating the differential equations 
could best be estimated by integrating the heat conduc- 
tion Eq. (5) numerically and comparing the results with 
curve (c) of Fig. 16. Again here the errors are probably 
not more than a few percent. The other physical 
limitations of the present model do not warrant a care- 
ful check of this point. 

We now pass to a discussion of some result of calcula- 
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tions determining observable effects of hypersonic trails 
of blunt bodies. 


(3.4) Thermal Radiation and Ionization in Equilibrium, 
Conduction-Controlled Laminar Trails 

The radiation intensity as a function of axial position 
along the trail at an altitude of 100,000 ft and for ve- 
locities ranging from 15,000 to 35,000 ft/sec is depicted 
in Fig. 26; the dashed curves give the intensity in the 
optical region, the solid that of the infrared. It should 
be noted that in order to get the absolute magnitude 
of the intensity as a function of the axial physical dis- 
tance from the p > p., points, in centimeters, both co- 
ordinates have to be multiplied by the nose radius (in 
centimeters) squared. Thus, for a given detectable 
threshold of radiation, large bodies may leave trails that 
last for miles. The radiation intensity in the optical 
region is a very steep function of flight velocity: it 
varies as the fourth power of velocity in the low-velocity 
range and as the tenth power in the high range. In 
the infrared range the previously mentioned velocity 
scaling is also applicable if the region close to the 
pb > p. point is disregarded for speeds higher than 
20,000 ft/sec. 
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From Fig. 26 it can be concluded that in interpreting 
experimental radiation data from the trail, it will be 
very important to determine the flight velocity with 
accuracy. Fig. 26 also shows that at low velocities 
(15,000 and 20,000 ft/sec) the energy per unit length 
in the infrared range is between one and two orders of 
magnitude larger than in the optical range. At higher 
velocities the trend is reversed. The general steep de- 
pendence of intensity on axial position should be noted. 
At velocities larger than 25,000 ft/sec, as the trail cools 
the NO radiation first increases and then decreases. 
The reason for this behavior is that at first the tem- 
perature in the trail is sufficiently high so that nitric oxide 


dissociates. As the trail cools NO forms up to a peak 
concentration and then decays. 

Fig. 27 gives the radiation intensity per unit length 
for an altitude of 250,000 ft and a velocity of 25,000 
ft/sec, from which it can be observed that the lower 
gas density at high altitudes drastically reduces the 
emitted radiation and increases the cooling rate. The 
trail will probably appear quite short for conventional 
detectors; thus, if it were of interest to estimate the 
optical radiation that decays drastically in about 5r,, it 
would be necessary to calculate the conduction and 
expansion cooling simultaneously when dealing with 
bodies up to about 30-cm nose radius. Under these 
circumstances the conduction-controlled trail is equal 
to or shorter than the expansion trail. Fig. 27 indicates 
that the optical radiation radically decays within a 
length equivalent to 5r, radii from the p —> p.. point. 
Thus for bodies with up to 30-cm nose radii the decay 
length would be 150 radii or less. Now it was seen in 
Section 2.4 that this is the length it takes for the pressure 
to decay from the value at the body to the free-stream 
value. Thus, the present result of Fig. 27 and the ex- 
pansion-controlled estimate (which could be obtained 
by the method described in Section 2.4) would be in 
error. Again, in order to obtain the correct estimate 
of the optical radiation, it would be necessary to include 
the conduction cooling in the expansion calculation. 
However, in the case of the infrared radiation (Fig. 27), 
this decays slowly and an error would be introduced 
only in considering bodies with nose radii of less than 
3 cm. Altitudes and velocities other than those 
given in Figs. 26 and 27 can readily be computed. We 
have only computed one condition at the higher altitude 
in order to investigate the usefulness of expressing the 
results in terms of the dimensionless variable X¥ = x*/ Ré 
introduced in Sections 3.2 and 3.3. This independent 
variable would be very convenient if it were not for the 
fact that, for example, at constant flight velocity a 
change in altitude affects the absolute magnitude of 
the thermodynamic variables in the trail. In order to 
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make the point, we have plotted in Fig. 28 the radia- 
tion intensity distribution (divided by a suitable power 
of the free-stream density) as a function of X for a flight 
velocity of 25,000 ft/sec at two different altitudes, 
100,000 and 250,000 ft. It can be seen that the axial 
dependence of the infrared radiation is correctly corre- 
lated by the parameter X. The altitude dependence 
of that radiation is well correlated by the 3/2 power of 
the density. However, in the case of the optical radia- 
tion, the parameter X does not correlate the results 
very well in the axial direction. The values of X for 
the optical radiation at 100,000 ft would have to be 
stretched by a factor of about three in order to match 
the 250,000-ft radiation distribution. The density de- 
pendence of the optical radiation seems to be well 
correlated by the cube power of the density. The cube 
power results from the fact that in the altitude-velocity 
range considered, the radiation is collision-limited.'! 
The total integrated radiation can be obtained by inte- 
grating Eq. (4) with respect to the axial coordinate, x—, 
ie., J, = f{ Idx. Since J,/r, ? is a unique function of 
x/r,”, it will be convenient to express this integral as 


I. (watts) ” I, (watts/cm) x 
= d -lem-* (20) 
0 Vn” 


7,* (em*) r,2 (cm?) 


n 


Fig. 29 was obtained by using Eq. (20) in connection 
with the results of Fig. 26. From Figs. 13 and 29, it is 
now possible to estimate the total thermal radiation 
emitted by the trail of a blunt body. The total radi: 
tion in the expansion-controlled trail or from the nose 
cap scales as 7,*, while the radiation of the conduction- 
controlled trail scales as 7,*.. The reason for this is that 
they both scale as the cross section of the trail, which 
gives for both an 7,” dependence; in addition, the 
length of the trail and the nose cap scale as r, in the case 
of the expansion-controlled and 7,” in the case of the 
conduction-controlled trail. Depending on the absolute 
value of the nose radius 7,, (in centimeters), either the 
expansion-controlled trail or the conduction-controlled 
trail will emit more thermal radiation. Fig. 30 shows 
the effects of nose radius (1 to 10? cm) on the total 
emitted thermal radiation; the results of Fig. 13 were 
included. 

Fig. 31 gives the electron density distribution on the 
axis of the trail for several velocities at an altitude of 
100,000 ft, and Fig. 32 gives a complete map of the 
electron density field for several flight velocities. The 
lowest electron concentration plot was for n, = 10° 
electrons/cm*, which corresponds to the typical elec- 
tron densities occurring in the ionosphere—electro- 
magnetic reflection and scattering from the trail by 
electron densities of this order of magnitude are con- 
fused with the ionospheric background. Again here, 
as in the case of thermal radiation, the electron density 
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has a very steep velocity dependence. The importance 
of accurately determining flight velocity in experimental 
investigations cannot be overemphasized. The length 
of the conduction-controlled electron trail (Fig. 32) is 
proportional to the nose radius squared and thus is quite 
large for blunt bodies. 

Fig. 33 gives the electron density along the axis at an 
altitude of 250,000 ft. At these altitudes the trail 
cools very rapidly, and ambipolar diffusion becomes 
important if the body is small. This effect will not be 
considered here. 

Fig. 34 is an attempt (similar to Fig. 28) to correlate 
the results of electron density distribution (non- 
dimensionalized with respect to the proper power of 
density) as a function of the axial variable X. The 
variable X seems to correlate the electron density dis- 
tribution, at a given velocity, for different altitudes 
within a factor of 10. 


(4) Concluding Remarks 


We have treated here a fairly restricted problem 
which may not be too realistic in technically important 
situations, the most important restrictions being the as- 
sumption of equilibrium and laminar flow and the fact 
that the existence of the boundary layer in the wake 
core has been disregarded. However, the blunt-body 
assumption is quite realistic. Nevertheless, the calcu- 
lations have led to an evaluation of the scaling laws in- 
volved in the problem. For high velocities and inter- 
mediate altitudes, the values depending on body size, the 
present calculation in the case of laminar trails gives a 
lower limit for the electron density distribution. 

The flight conditions for which laminar-turbulent 
transition will occur in the trail, as well as its rate of 
growth, will have to be determined experimentally. 
When this is done the present calculations will serve as 
background for extending the analysis to the turbulent 
trail. 


Appendix 


Calculation of the Length of the Frozen Wake 
By assuming that a single diatomic gas approximates 
the behavior of air (which is done only for order-of-mag- 
nitude estimates) and using typical chemical recom- 
bination rates k, for the chemical species in air, it is 
possible to calculate the length of the trail that will be 
chemically frozen. For the three-body recombination 
of two atoms of species A by means of a catalyst M 
leading to a molecule A, and an energetic catalyst M*, 
the reaction can be written as 
kr 


A+A+M->A,; + M* 
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where k, is the recombination rate coefficient. The 
kinetic differential equation is given by 
d{A] 


— 2k,[A]?[M] 

dt | 

where the brackets denote concentrations in gram- 
moles per cm*. Using a value of k, of 

4500 
T (°K) 


k, = 10" | cm® mole~? sec~! (A.1) 


one can estimate a recombination time ¢, from 
t, = }2k,[A][M]}-! (A.2) 
The flow time ¢; can be approximated by 
tr = 1;/ta (A.3) 
where /, is the length of the frozen wake and w., is the 
flight velocity. For frozen flow, we will assume that 


the flow time is small compared to the relaxation time. 
Thus, for our purpose let 
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Assuming equilibrium flow in the forward portion of the 
blunt body, the temperature and chemical composition 
can be determined. Then Eqs. (A.1) through (A.4) 
will permit the evaluation of /,; for any flight velocity 


and altitude. 
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Nonsimilar Solutions of the 
Compressible-Laminar Boundary-Layer 
Equations With Applications to the 
Upstream-Transpiration Cooling Problem! 


ADRIAN PALLONE* 
Aveo Research and Advanced Development Dwision 


Summary 


\ new method is presented for predicting the boundary-layer 
characteristics downstream of the porous region of an injection- 
cooled surface. The method consists of a general scheme for 
obtaining nonsimilar the 
boundary-layer equations and is formulated along the following 
The viscous domain is divided into N curvilinear strips 


solutions of compressible-laminar- 


lines. 
The governing equations are then integrated along the coordi- 
nate normal to the body from the surface to the boundary of each 
As a result, one obtains a set of independent integro- 
by ex- 


strip. 
differential relations. 
pressing the integrands as polynomials, the coefficients of which 
are functions of the unknown values of the velocity and tem- 


The integration is carried out 


perature on the strip boundaries as well as of the imposed bound- 
ary condition at the wall and at the outer edge. After the 
integration is performed, a set of ordinary first-order differential 
equations is obtained. The set of equations may be solved for 
given initial conditions by a numerical integration scheme such 
as the Runge-Kutta method. Several numerical examples of 
interest are presented. 


Symbols 


. . . ~ . Tr 
( = local skin-friction coefficient = 
( pit,* 2) 


= ( pyuyx/pyc)~/2 2a,(t/d)"!2(1/c) 
p ad g 


ce = specific heat at constant pressure 

f(0) = variable proportional to the blowing rate in 
reference 11 

Fix, , Fs, = integrals defined after Eq. (12) 

G = function of 4, m4, and wy, cf. Eq. (24) 

H = stagnation enthalpy = (u?/2) + c¢,T 

l = nondimensional ratio pu/pimi 

L = reference length 

M = Mach number 


N = number of strips into which viscous region is 
divided 


p = static pressure 

Pr = Prandtl number c,u K, where K = con- 
ductivity 

O = Surface heat-transfer rate 

R = gas constant 
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R = gasconstant of mixture 

Re, = Reynolds number based on properties at 
outer edge of the boundary layer = 
UX pi / Mi 

Ss = Sutherland constant, 216°R for air 

t = transformed variable, cf. Eq. (9) 

T = absolute temperature 

T = an average temperature = (7). + 7,)/2 

u = velocity component in x direction 

i = velocity component in y direction 

x = coordinate along the surface 

y = coordinate normal to the surface 

7 = ratio of specific heats 

n = transformed variable, cf. Eq. (10) 

6 = boundary-layer thickness 

A = nondimensional boundary-layer thickness = 
( pitt; L/ my )(6 2 /L?) 

Mm = coefficient of viscosity 

‘ = nondimensional space coordinate = x/L 

p = mass density 

T 1(0/On)(u/u;) 

Subscripts 

a = adiabatic wall 

L = value at the end of porous region 

r = reference value 

s = stagnation value 

t = value in the x-t or &n plane 

Ww = value at the surface (y = 0) for any value of x 

1 = value at the outer edge of the boundary 
layer 


= value at free stream 


(I) Introduction 


A“ EFFICIENT METHOD of protecting a surface from 
the high-energy gas stream consists of artificially 
thickening the boundary layer and altering its profile 
by the mechanical injection of gas normal to the solid 
surface. Flight vehicle design considerations usually 
dictate that the mass be added over discrete regions, 
since continuous mass addition involves the use of a 
porous surface material, which has reduced strength 
properties. 

The amount of coolant introduced by the localized 
mass-addition scheme is considerably greater than is 
necessary to protect the area of injection itself; the 
excess coolant is utilized to protect some downstream 


area. The persistence, downstream from the point of 
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injection, of the alteration of the boundary-layer 
characteristics determines the effectiveness of the cool- 
ing scheme based on such discrete mass addition. 

An analytic evaluation of the insulating effect of gas 
injection into the laminar boundary layer was first 
presented in reference 1. The analysis is based on the 
Karman-Pohlhausen integral scheme; however, addi- 
tional parameters were introduced in the profile which 
were determined by making the flow of mass momen- 
tum and energy continuous at the interface between 
the porous and nonporous sections. By handling the 
discontinuity of injection in this manner the profiles 
are made too “‘rigid,’’ since no decay mechanism is pro- 
vided for the additional parameter. A similar investi- 
gation using the Karmaén-Pohlhausen approach was 
presented in reference 2, where the additional parameter 
and continuity conditions of reference 1 were replaced 
by a continuity in wall shear. 

In reference 3 a numerical solution of the same prob- 
lem was obtained for the case of zero pressure gradient. 
This solution is presumably correct since no arbitrary 
assumption was made and the exact partial differential 
equations were integrated by a finite difference scheme. 
Results obtained by the three methods, for a particular 
case, were compared in reference 3. The finite differ- 
ence solution lies between those of reference 1 and 2, 
and there are large discrepancies between the results 
of the three methods. It can be concluded, therefore, 
that except for a very rough qualitative estimate of the 
effect of the upstream injection on the boundary layer 
characteristics, a solution for this problem based on 
the Karman-Pohlhausen approach is unsatisfactory. 

In this report a new method of obtaining nonsimilar 
solutions of the boundary-layer equation with and 
without a streamwise pressure gradient is developed 
and applied to the problem at hand. The method 
combines the Dorodnitzin integration scheme‘ with the 
Pohlhausen approach and is formulated along the 
following lines: the viscous domain, from the body 
surface to the outer edge of the boundary layer, is 
divided into N curvilinear strips. The boundary-layer 
equations are then integrated along the coordinate nor- 
mal to the body from the surface to the boundary of 
each strip. In order to carry out the integration, the 
profiles are represented by polynomials; in contrast to 
the Karman-Pohlhausen method, however, the co- 
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efficients of these polynomials are not known a priori 
but are expressed as functions of the unknown values 
of the integrand on each of the strip boundaries. 

After the integration is performed a set of ordinary 
first-order differential equations is obtained. The 
values of the unknown function at the boundary of each 
strip may be obtained by integrating the set of equa- 
tions along the streamwise direction from given initial 
conditions. The known function may subsequently 
be used to determine any other desired boundary-layer 
characteristics. 

The method of solution is developed for a dissociating 
gas in thermodynamic equilibrium; however, the illus- 
trative examples in this paper pertain to an ideal gas 
only. Numerical results for a dissociating gas will 
be reported in future work. 


(II) Basic Equations 


Eqs. (1) through (4), describe the two-dimensional 
boundary-layer flow of dissociated air in thermodynamic 
equilibrium with a Lewis number of unity and a con- 
stant Prandtl number. 

Momentum equation in x direction: 


puu, + piu, = —pr + (uty), (1) 
Continuity equation: 
(pu), + (pv) y = 0 (2) 


Energy equation: 


pul, + po, = 


(F%) +[e(-s)S)] © 


Equation of state: 
pRT (4) 


Egs. (1) through (4) give p, u, v, and H as functions of x 
and y for a given set of boundary and initial conditions. 
The quantities at the outer edge of the boundary layer 
are assumed known from inviscid flow consideration. 
A convenient method of solving the above-formulated 
equations is to combine the numerical integration 
scheme of Dorodnitzin‘ with the Pohlhausen approach. 

Egs. (1) and (3) are multiplied by dy and integrated 


from y = 0 toy = »,, wherek = 1,...,N. Thus 
one obtains: 

yk pu dp ( =) ( ~) ‘ 
7 ase dy = a = a co (5) 
pitty dx Jo 7 OV] w Ov]; 
u 

Pf dy = 

Pilly 
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The following boundary conditions are applicable to the equations: 
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at y = OU: u=v = J, 


at y = 6(x): 


H = H,, or 0H/Oy = const. (7) 
u = uy, Ou/Oy = 0, H = NM, 0H/dy = 0 (8) 


Eqs. (5) and (6) can be written in a more convenient form if the Howarth-Dorodnitzin transformation is introduced : 


t = 


Moreover, 


[VN — (k — 1)]/N. 


where 7, is chosen as 7, 


t, =| 
0 


Thus one obtains: 


- dy (9) 
Pl 


: dy = nb: (10) 


(Fyd'/2) + AMF! + Fig} (In (pret) 2)’ + (in wy)! — AFoy(uy/ttr)’ + AFP (Inu)! — te + 7 0 (11) 


and 


(F)’ 2) + AF 4,’ + NF, [In (py2¢1p3)! 2]’ —_ A Fox (Hi, / FT;) 


i fi l (HH l u,* O fu \* 
— _ (Pr — 1) (12) 
Pr Fi, nu Pr Hi, nk Pr 2H, On uy k 


where: 


‘ nk f Uy u* 
Ky, = -_— — = dn 
0 Uy Uy 4° 
‘ f py Uy, U 
( —- dn 
p Uy Uy 


r= {" u (F ye 
Tk a a 


)’ denotes differentiation with respect 


fT 
aa 


and where ( 
to § =<x«/L. 

Polynomials in are assumed for u/m, and H/H;, and 
the coefficients of these polynomials are selected so that 
appropriate boundary conditions are satisfied, while 
the remaining unknowns are determined so that Eqs. 
(11) and (12) are satisfied. The essential difference 
between the Karman-Pohlhausen method and the 
method used here is that in the former the partial dif- 
ferential equations are satisfied only in the average over 
the entire boundary-layer thickness—i.e., N = 1, or 
k = 1 only—while here the partial differential equa- 
tions are satisfied in the average also over additional 
subintervals of the boundary-layer thickness (NV > 1). 

Additional accuracy may be realized if one follows 
the scheme customary in refining integral methods—.e., 
if one assumes the polynomials to be of sufficiently high 
order that they satisfy some additional condition at the 
wall and at the outer edge of the boundary layer corre- 
sponding to those which an exact solution of Eqs. (1) 
through (4) would automatically satisfy (cf. reference 5). 

We assumed polynomials of degrees N + 2 for u/m, 
and of degree NV + 3 for H/A. 


N+2 N+3 
u N+2 H ; 

-_ = > ann”, _ + & bn” ( l 3) 
uy m=0 HN, m=0 


The additional coefficient in the thermal profile is as- 


. y¥-—1 nm f HT, =) 
Fy. = | 1 - M,? = 1 
( T 2 yf & uy" w 


pylt; Lé?? 


Mi i 


sumed to serve as an unknown in addition to the com- 
mon boundary layer thickness, and is an alternative, 
for example, to a separate thermal boundary-layer 
thickness (cf. reference 5). 

The boundary condition in the £-» plane correspond- 
ing to those in the x-y plane given by Eqs. (7) and (8) 


aves 
’ oe , Il H, (“") ) 
a = (0): = (), : ) rive 
. ui, x & Gi” 
(14) 
u H 
aty = 1: =1=—>=1 | 
uy i, } 


The additional boundary conditions can be obtained 
from Eqs. (1) and (3) at 7 = 0: 


non - 2[12(t) 
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re) O/H — f0u/u,\? 
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By means of Eqs. (13) through (16) expressions for the F integral can be obtained, and Eqs. (11) and (12) can be 
rewritten as a set of 2/N first-order differential equations. 

In order to illustrate the method of solution, we will discuss in detail the case of an ideal gas; the procedure for a 
gas in thermodynamic equilibrium will be outlined in the Appendix. 


(A) Ideal Gas 
This case is further simplified by assuming a linear viscosity-temperature relation in the form introduced by 
Chapman.® Thus: 


l= e(Tur)s/(T re); C= ul y,/Tux, - (T., Ty)? (Tz, + S)/(T + S) | (17) 


The F integrals for the ideal gas case become: 


. Uy. V+2 Om m+1 N+2 Qin? 2m +1 ” N+2 inne m+2 
Fy. = - — a + 2  B 9 + 
U1 may M+ 1 mal om+ 1 ma2 mM+2 
9 A? Gm ** yA? GnwH18mm OV? 
20 9 +... +2 ps 2(N + 2) 
m=3 m = 0 m=N+2 2\i T 4 
_ re 
Foy = Ah? OmNk ail 
m=1 ™ -~ 7 


¥ ‘ = l Be Din N+2 Om wane 
F 3 = F _— F =_ ( l + = M 2 m +t 7 + 
, , . . 2 2», m = 1 ™ > 2m + ] 


) ~ QAmNx salutes 9 N+2 Q2OmNk wre N+2 Q(nt+ YD OmNk 2(N +2 
22> — +2 > St DO — Fy 
m=2 m+ 2 m=3 mM+3 m=N+2 2(N + 2) 
H, N+2 OmNk m +1 N+2 OmNk m+1 V+2 OmNk m+2 
Fu = — (oS ee 45 
HAH, 2, m+) mal M+ 1 mol M+ 2 
N+2 m +3 N+2 m+N +4 
- amNk “ Ams 
b > +...+tbyi3 > 
m=1 ™ + 3 m=1 ™ + N + 4 
wherekR = 1,...,WN. 


The coefficients a,, and b,, can be obtained in terms of u,/u; and H,/H, respectively by solving the following set 
of algebraic equations: 


V+2 AD N+3 
Up k = 
yon » AmNk”; H _ ss Drnx” k = | ae cs N 

1 m=0 1 m=0 


subject to the conditions in Eqs. (14), (15), and (16). 

The solution of the system of Eqs. (11) and (12) for the zero pressure gradient case will be discussed first. This 
solution is of great value in evaluating the method since it can be compared with existing solutions, in particular 
the exact solutions of reference 3. 


(1) Zero Pressure Gradient—For the flat plate case, Eqs. (11) and (12) become: 


tidy a wn O@fm N+2 . 
(; dé ais 5) ra ae dé (“*) (« a > — )- ' — 
1 dx —a@ : os. 1G eee 1 N+3 ” 
(; * tie 2 ellie iy ~ Pr E “2 J ‘ 


1 — Pr 2M,-2\—! u, N+2 . 
( Ne a ) 2 < > man; | 0 (19) 
Pr 7 ] U1 m=] 
where \ = A/c. 


In order to complete the system a set of initial conditions and the number of strips V must be specified. Since the 
accuracy of the result will depend on the number of strips prescribed, the solution will be carried out for several values 
of N. 

Eqs. (18) and (19), with the accompanying initial and boundary conditions can now be solved for the 2N un- 
knowns, namely, u;/2, \, H,/H, (k = 2, 3,..., N) and either 77,,/H, or [(0/0n)(H/Hj) |, depending on whether 
the wall temperature or the rate of heat transfer is specified. 

It is interesting to note that for the particular case of Pr = 1 and for a given constant H/,, equal to H/,, of the porous 
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region, the energy Eq. (3) is satisfied by the Crocco integral. 
velocity field by the following expressions (cf. reference 7) : 


H = (MW, — H,)(u/m) + He 


Flow With Pressure Gradients 
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The thermal field is thus directly related to the 


Q = (77; = H,,)(pittx? 2)cy (20), (21 
In this case, Eqs. (11) and (12) can be written in the following form: 
ge l N+2 
+A Fy —Inw — a— p M Onn.” 0 (22) 
dé c m l 


LU 


~ mbar) — 
2M N+2 
(1 — ab : 5 yy (2! 2 Die MA mn,” ) =(Q (23) 


Ui m=1 


where G = (pitty)! 2 (24) 


Eqs. (22) and (23) constitute two coupled systems of equations to be solved for the unknowns x;,, \, J7, and either 


i HT, or [(0/0n)(H/Aj) |. 


(B) Numerical Examples 

In this section, results of the above analysis are illus- 
trated by several examples. 

The general problem consists of determining the 
boundary-layer characteristics downstream of a par- 
tially injection-cooled surface immersed in an ideal gas 
stream, as shown in Fig. 1. In the solution, the ve- 
locity and thermal profile as well as the boundary-layer 
thickness and the wall temperature were assumed 
continuous at the junction of the porous and nonporous 
region. 

The numerical solutions were obtained by means of 
two programs. In the first program the governing 
equations were cast in terms of the dependent variables 
u,/u,, \, H,/H; and either H,,/H, or [(0/0n)(H/H,;) |, 
The accuracy of these equations was improved by per- 
forming the calculation in “‘double precision”’ and re- 
taining the “‘most significant half’ for the coefficients. 
Even though the machine running time increases 
roughly tenfold when double precision is used, the sets 
of equations for NV = 3, 4, and 6 were obtained in ap- 
proximately 20 minutes. 

In the second program, the set of first-order ordinary 
differential equations corresponding to the value of NV 
of interest was solved from given initial conditions by 
use of a Runge-Kutta fourth-order method. In this 
program, double precision was not required. The 
integration step At was allowed to vary from 0.025 to 
0.1 along the body according to an arbitrary Ac, cri- 
terion. The running time for the pressure gradient 
case N = 6 for a total Ag of 4 was approximately 20 
minutes. This running time could be considerably re- 
duced if the additional boundary condition (15) were 
relaxed, which would permit uncoupling of the two sets 
of equations at each station. The accuracy of the re- 
sult with condition (15) removed must, of course, be 
checked. 

(1) Zero Pressure Gradient 
for zero pressure gradient the surface is considered flat. 
Coolant air is assumed to be injected normally to the 
surface at a nonuniform rate such that the wall tem- 
perature at a distance downstream of the leading edge 


As a numerical example 


Fic. 1. Schematic of upstream injection technique 
is maintained constant. The remaining portion of the 
surface is assumed insulated. The following free- 
stream conditions are considered: 

M. = 3) T.. = 390°F Tw, = 0.57, 
M.=6) Pr=0.72 p. = 472 psf 
The skin friction and the recovery temperature were 
(18) and (19) by the 
velocity and thermal 


obtained by integrating Eqs. 
Runge-Kutta method from the 
profiles immediately downstream of the injection re- 
gion. 

The integration was carried out for values of V equal 
to 3, 4, and 6, and for the set of initial conditions ob- 
tained from reference 11 (see Table 1). Figs. 2 and 3 
show the ratio of the local skin-friction coefficients for a 
plate with and without injection but with identical wall 
temperature distributions. Also shown in Fig. 2, for 
purpose of comparison, are the results in Fig. 7 of refer- 
ence 3. It is noted that excellent agreement is ob- 

TABLE 1 


Initial Profiles™ 


4 HH + Hy / 4/e 
‘ ios) (Wee  r | «| |See 
ir 1.¢ 1, 000 1. 00 1 0 1. 00 
9937 1. 000 00 oc 8 98¢ 1.0 
0. 9807 0. 9887 931 96 4 0. 8841 4 0.938 
4 8056 0. 884( 0. 8863 0. 7354 0.834 8 7428 
( 6959 4358 64 4 284¢ 616 ¢ 4 
4i¢ 0. 565 4 1 8 4784 4) 0.46 0.429 
‘ 
0 ¢ 00 0 
81 9964 9891 0. 9956 9565 0.9821 0. 9884 
8 0.884 886 834 829¢ 8 0.7428 c 
408 0.1718 1.5039 0. 4686 
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Fic. 2. Ratio of the local skin-friction coefficients for a plate 
with and without injection, but with identical wall-temperature 
distributions; f(0) = —1.0. 








Fic. 3. Ratio of the local skin-friction coefficients for a plate 
with and without injection, but with identical wall-temperature 
distributions; f(0) = —1.0, —0.75, —0.50. 








3 4 


Fic. 4a. Comparison of skin-friction results for N = 3, 4, and 6; 
f(0) = —0.5. 








Fic. 4b. Comparison of skin-friction results for N = 4 and 6; 
f0)= —1.0. 


tained between the present results and those of the 
finite difference solution. The discrepancy between 
the present result and those of references 1 and 2 is 
due primarily to the inherent weakness of the Karman- 
Pohlhausen method. The inadequacy of this method 
is manifested by the fact that only one initial condition 


1961 


from each initial profile is required, which greatly 
weakens the influence of the initial profile on the down- 
stream development of the boundary layer. The 
authors of reference 1 tried to remedy this situation by 
introducing additional parameters into the downstream 
profiles, obtained by making the flux of mass momen- 
tum and energy continuous* at the interface between 
the porous and nonporous region. This technique 
made the profiles too “‘rigid,’’ since it did not provide 
any decay mechanism for the additional parameters. 
In reference 2 also an attempt was made to enhance the 
influence of the initial profile on the downstream solu- 
tion by making the velocity profile of sufficiently high 
order that an additional differential equation relating 
(0u/Oy),, 6and x was required. Although this method 
requires two initial conditions—namely, the bound- 
ary layer thickness and the wall shear at the injection 
region—the initial profile has still only a weak influence 
on the development of the solution downstream. 

In Figs. 4a and 4b the skin-friction results obtained 
by the present method for several values of NV are pre- 
sented.** For the small injection rates [f(0) = —0.5| 
the values obtained for NV = 4 are within 2 percent 

Errors on the order of 
For the larger 


of those obtained for NV = 6. 
6 percent are noted for the NV = 3 case. 
injection rates [f(0) = —1], the discrepancy between 
the NV = 4 case and the N = 6 case increases to 5 or 6 
percent, and for the NV = 3 case large errors are noted. 
(No curve is shown for N = 3 in Fig. 4b.) 

Figs. 5a and 5b show the variation of surface tem- 
perature downstream of the injection region of a flat 
plate at M/.. = 3 and 6, respectively, for two values of 
f(0)(—1 and —0.5); N was assumed equal to 6. For 
comparison the result of reference 3 obtained by finite 
difference solution is also shown in Fig. 5a. A sub- 
stantial reduction in recovery temperature is realized. 

(2) Favorable Pressure Gradients—As a_ second 
example, quantitative effects of upstream injection are 
obtained for the curved surface shown in Fig. 6; a nega- 
tive pressure gradient is induced. 

The leading 20 percent of the surface is assumed to 
be flat and porous. Air is assumed to be injected at a 
nonuniform rate so as to give a constant wall tempera- 


ture. The following free-stream conditions are con- 
sidered: 
M. = 6, Ty, = 0.7T7,, T. = 390° R 


Pr = 0.72, P. = 472 psf 


The inviscid flow quantities such as M4, m and 7; 
were determined by the shock-expansion method 
neglecting reflected waves. The boundary layer char- 


acteristics were obtained by integrating Eqs. (22) and 


* One may note that this condition is automatically satisfied 
in the present method since the profiles are assumed continuous 
at the interface between the porous and nonporous section. 

** The discrepancies for the various values of N are due pri- 
marily to the fact that the initial profiles are not exactly repre 
sented for the small values of N. 
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TABLE 2* 


Initial Profiles!! 


* Mt \ 

006 1. 0000 1. 000C 1. 0006 

7 OC l 1. 0000 ( 87 o0oc 
0.9517 0. 9876 0.9315 0.9732 0. 8841 77 

0. 805¢ 0. 886 0. 735( 0.8298 5983 0. 7265 

71 0.698 0. 4358 0.627 0. 2846 24 

0. 241¢ 0.54 72 48.28 ( 4 434 

I | | | | ; 
*M 5 2 is the Mach number obtained behind the leading-edge shock 


corresponding to a free-stream Mach number of 6 


(23) for N = 6 and for the sets of initial conditions 
obtained from reference 11 shown in Table 2. 

In Fig. 7 the values of the local skin-friction coeffi- 
cients along the favorable pressure gradient surface are 
presented. The coefficients are given in terms of the 
coefficient of the local flat plate without injection but 
with the same wall temperature distribution as for the 
curved surface. For purposes of comparison a curve 
showing the variation of skin-friction coefficient for the 
case f(0) = —1.0, obtained by neglecting the pressure 
gradient (local flat plate), is also included. From this 
comparison one may note the strong effect of the pres- 
sure gradient on the rate of skin-friction recovery. 

The recovery wall temperature for the favorable 
pressure gradient surface is shown in Fig. 8. Again for 
purposes of comparison a curve obtained by neglecting 
the pressure gradient is included. One may note 
that the pressure gradient affects the recovery tem- 
perature only slightly, but in the direction opposite to 
the one obtained for the skin friction. A similar result 
was also shown qualitatively in reference 1. 

Figs. 9a and 9b show a typical development of the 
velocity and stagnation enthalpy profiles downstream 
of the injection region for the case f(0) = —1.0. The 
value of \ for the same case is shown in Fig. 10. From 
this figure it is interesting to note that as the velocity 
profile becomes fuller, the value of X decreases to a 
minimum at a short distance downstream of the injec- 
tion region and then starts to increase at a rate which 
becomes almost linear with &. 


(III) Conclusions 


The boundary-layer characteristics for flows with 
and without pressure gradient can be determined for 
any given initial conditions by integrating numerically 
the set of first-order ordinary differential equations 
developed here, provided the Lewis number is unity 
and the Prandtl number is a given constant. 

The method is applied to the calculation of the local 
skin-friction coefficient and wall temperature dis- 


Fic. 5a. Surface-temperature distribution for an adiabatic 
wall downstream of a transpiration-cooled region; M, = 3, Tw, 
= 0.5 7.. 








4 


Fic. 5b. Surface-temperature distribution for an adiabatic 
wall downstream of a transpiration-cooled region; M, = 6, Tx 
= 0.5 T.. 






POROUS 
REGION 


0.02 -— 


Fic. 6. Curved surface inducing a favorable pressure gradient. 





4 8 22 26 


E 
Fic. 7. Skin-friction coefficient for a surface inducing a favor- 
able pressure gradient in terms of a local flat plate without injec- 
tion, but with an identical wall-temperature distribution; 
Moa = 6, N = 6. 
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pi/p = (pi/p,)[0.8651(h/h,) — 0.0053(h/h,)?] (A-1) | 


l = (p,/u,/ pis) [1.5852 (h/h,)-/? — 
























(0.2829 (h/h,)-!]  (A-2) : 
where h, = 118 Btu/lb, nw, = 1.7235 X 10-4 g : 
oat | cm/sec, pp = P\/RT,, R = 0.6886 XK 10° cal./g | 
=-= K, 7, = 273.16 °K, and 
ot oe Coa h H,H 1 (u2\/u\? 
Fic. 8. Surface-temperature distribution for an adiabatic curved h, o E H, 7 5 (“! ye) | | 
surface (favorable pressure gradient); Mo = 6, N = 6. I 
The Prandtl number has been assumed constant | 
throughout the analysis; this assumption is justifiable ! 
for both ideal and real gas cases, in view of the result Sid 
shown in reference 10. The value of the constant may 7 ; 
range from (0.7 to 0.82 depending on the problem at si 
hand. an¢ 
The relations for F,, /:, and F; are the same as those ' 
given for the ideal gas case. However, the expression the 
for F; must now be obtained by substituting Eq. (A-1) reg 
. into the F; integral expression. 
(Continued on page 492) or 
in 
mu 
Thi 
fe) tral 
tur 
Re} 
con 
dep 
I 
u/u, ree 
Fic. 9a. Comparison of velocity profiles for various stations aks 
downstream of the injection region (pressure-gradient case); ; 
(0) = -1.0. ~ 
ol 
lay 
tribution downstream of the injection region of a par- ™ 
tially porous flat plate. my 
The numerical results, obtained for the particular eas 
case where the surface downstream of the injection flow 
region is flat, have been compared with the exact re- aes 
sults of reference 3 (obtained by utilizing finite differ- — 
ence techniques). Complete agreement is noted when H/H, sa 
N is taken equal to 6. The skin friction coefficients Fic. 9b. Comparison of stagnation enthalpy profiles for vari- i, 
for N = 6 and N = 4 agree within 2 percent for the ea age haat" all Se ee ae a 
small injection rates and within 6 percent for the large ” 
injection rates. . F the 
' ste 
Appendix A a 
Herein is presented an outline indicating the pro- ue 
cedure to be used for a dissociating gas in thermody- = 
namic equilibrium.* The physical properties of the 
gas, appearing in Eqs. (11) and 12) as Pr, /, and p/pu, 
can be specified as functions of H/H, and u/im, by I 
following the procedure of references 8 and 9. Suffi- 196 
cient accuracy may be realized, in the range of h/R7T) < 
220 and.of P; from 0.1 to 10 atm., by the following is 
analytical expressions: pee 
Se sea a ry ey ey , § 
* Numerical results for this case will be presented in future Fic. 10. Nondimensional boundary-layer thickness parameter on 
work. versus & (pressure-gradient case). § 
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An Investigation of Separated Flows 
The Pressure Field 


Part I: 





A. F. CHARWAT,* J. N. ROOS,** F. C. DEWEY, Jr.,§ anv J. A. HITZS§$§ 
University of California at Los Angeles 


Summary 


The present article describes an investigation of several types of 
separated regions such as blunt-base wakes and cavities formed 
in cutouts in the boundaries and ahead of or behind two-di- 
mensional steps in supersonic (Mach numbers 2 to 4) and sub- 
The for the 


and the pressure field are described in this paper 


sonic flow conditions the existence, geometry, 


\ second article (to be published ) will describe investigations of 
the internal flow and the heat transfer across such separated 


regions, 


It is found that there is a maximum (critical) ratio of the length 
of the separated free-shear layer to the depth of the depression 
leaving 


in the boundary beyond which the cavity collapses, 


mutually independent separated regions at each protrusion 
This critical length changes greatly upon laminar-turbulent 
transition in the oncoming boundary layer; in either laminar or 


Mach 


A semiempirical correlation predicting the 


turbulent flow it is approximately independent of and 


Reynolds numbers 
conditions under which the flow will span a depression of arbitrary 


depth is proposed 


Detailed distributions along the boundaries of a 


cavity (in turbulent flow) are presented as a function of the ratio 


pressure 


which is found to be 
the pertinent similarity parameter. notches (L/L, 
<0).5) the impact pressure due to the reversal of the inner portion 


of the cavity length to the critical length, 
For short 
of the shear layer at recompression tends to thicken the shear 


layer and a type of boundary layer-free stream interaction 


governs the pressure field. The pressure in the cavity is nearly 
constant and can be higher than free-stream. In long notches 


UF 0.5) 


curves back gradually ahead of the recompression point. 


the shear layer bends inward at separation and 
The 
floor-pressure variation is pronounced and the recovery pressure 
at reattachment is small. The variation of the drag coefficient 
with Mach number reflects the change from one to the other 


mechanism of recompression 


surveys of the Mach-number distributions in a 


wake and the mixing region behind its throat, as 


Detailed 
blunt-body 
well as in the shear layer spanning a cutout in a wall, are presented 
and analyzed. It is found that, in general, the assumptions of 
the simple supersonic-wake models which rely on a principle of 
steady flow with mass conservation in the cavity are not adequate 


for cavities in which there is recompression against a boundary 


Results showing the influence of the thickness of the initial 
boundary layer (in the range of 0.3 to 3 times the notch depth) 
and of the geometry of the notch are also presented. 
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Introduction 
v I ‘ue stupy of flows containing wakes and separated 
cavities is one of the oldest problems in aerody- 
There has recently been a growing interest 
This is 


namics. 
in the heat-transfer properties of such flows. 
intimately connected with the structure of the flow 
inside the cavity, about which almost nothing is known. 


The most complete conceptual model of these mixed 
flows currently available was first proposed by Crocco 
and Lees! and restated in a simpler form by Chapman 
et al.” and Korst et al. The problem is treated as an 
interaction between the external potential flow and 
the cavity proper. The profile of the viscous free- 
shear layer provides a link between these regions 

However, the flow within the cavity proper is out- 
side of the scope of the shear-layer models. There 
is ample evidence that this region is not “‘dead,” 
as postulated in the theory. It is a complex, vortical, 
probably three-dimensional flow and, as indicated by 
a growing body of experimental evidence, often an un- 
steady flow. As long as it is dynamically ‘“‘weak”’ in 
comparison with an external supersonic field, the mean 
geometrical characteristics of the wake and, hence, 
the mean pressure in it, are not strongly influenced by 
its internal structure. For this reason analyses which 
neglect it are successful in correlating base-pressure 
data. However, it is questionable whether this re- 
mains true in connection with transfer properties across 
the region, that is, whether the conduction (of heat, for 
instance) through the shear layer is the controlling 
step in the process or, indeed, an essential step. For 
instance, should there be a periodic unsteady mass 
discharge from the cavity, the conduction through the 
shear layer could turn out to be of secondary im- 
portance. 

The present research was undertaken for the purpose 
of studying the influence of as many parameters as 
possible on the internal and the external flow field in 
and about cavities. Several different types of separa- 
tions were used. Where necessary, observational tech- 
niques giving a clearer physical picture of the flow 
were stressed over quantitative results, and diversi- 
fication was chosen over the quantity of measurements. 

It was deemed useful to begin by a detailed study of 
the occurrence and overall structure of cavities, as a 
function of their geometries, the Mach number, the 
oncoming boundary layer and so on, which are dis- 


cussed in the present article. This is an extension of 
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Fic. 1. Pictures of cavity closure and typical pressure dis- 
tributions. M = 2.78, 6/H, = 0.44, turbulent flow, H,2 = H, = 
0.25 in. Sketch of models and definition of symbols. 
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Fic. 2. Critical closure distance for notches. 


some exhaustive measurements already in existence on 
simple separations behind and ahead of steps (for in- 
stance, reference 4). A second part of the study, to 
be published, will describe heat-transfer measure 
ments, as well as investigations of the steadiness, mag 
nitude, and direction of the flow inside the cavity. 


Equipment 


The work was carried out in two tunnels, one of which 
was instrumented for the measurement of the pressure 
field while the other was used principally for heat- 
transfer research. The major part of the data dis 
cussed in this paper was obtained in the first tunnel. 
This is a 3 X 3-in. test section, room stagnation tem- 
perature, variable stagnation pressure, continuous- 
flow unit equipped with profiled nozzle blocks for JA/ = 
2.78 and 1.9. The various test models are described 
in the text; most of the data was obtained on variable 
length and depth notches cut in the test-section wall 
(see Fig. 1b). The Reynolds number used in reporting 
the data is based on the measured boundary-layer 
profile just ahead of the notch (always turbulent), 
and the equation: 


6/S = 0.37 Re,-5 = 0.37 (UdS/v..)~/5 (1) 


The measured profile fits the 1/7-power velocity dis- 
tribution to which Eq. (1) applies. Large-aspect- 
ratio wedges, with and without splitter plates, mounted 
from the side walls (no sting) were also used. In these 
cases the Reynolds number is based on free-stream con- 
ditions and the wedge chord, but boundary-layer 
traverses were made to check (in turbulent flow) 
the approximate validity of Eq. (1). All pressure data 
were recorded on micromanometers to within better 
than 1 percent of free-stream pressure; all runs were 
repeated at least twice and considerable care was ex- 
ercised to avoid spurious effects such as leaks into the 
base region and other disturbances. 

Some additional pressure data (subsonic, JJ = 
and 3.5) were obtained in the heat-transfer tunnel. 
This tunnel is less thoroughly instrumented for pres- 
sure measurements (it reads to 5 percent of free-stream 
pressure at IJ = 3). It is equipped with a wedge- 
nozzle of small divergence so that there existed a 
static-pressure gradient over the test area (3.8 to 2.0 
percent of the static pressure at separation, per inch, 
in the range of Mach numbers tested). The pressure 
distributions obtained in this flow are reported in 
terms of ratios to the /Jocal free-stream pressure at the 
pertinent streamwise station. We believe that the 
pressure gradient was sufficiently small not to influence 
the flow field in the cavity. The Mach number and the 
Reynolds number Re, (based, as before, on the measured 
boundary-layer thickness) were evaluated at conditions 
corresponding to the separation point. 


») 


The Critical Length of Cavities 


A supersonic flow over a cut-out in a boundary (see 
Fig. 1) can exist in two distinct stable configurations 
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Fic. 3. Length-height ratios for separation wakes. (Re, based on undisturbed flow ahead of the separation.) Numbers in circles 
indicate references. Crosses are intercepts crossplotted from Figs. 2 and 6. 


When the length-to-height ratio of the cut-out is suf- 
ficiently large, the flow reattaches to the floor. There 
exists then a separated region behind the downstream- 
facing step (a separation wake) and another distinct 
and independent separation ahead of the upstream- 
facing step (a recompression wake). We term this 
a flow with a ‘‘closed’’ cavity. If the length of the 
notch decreases below a critical length, the flow bridges 
the cut-out entirely. This is called flow with an “‘open”’ 
cavity. 

Fig. 2 shows measurements of the critical closure 
length for various notches at several Mach numbers, 
differing by the shape of the recompression step and the 
height ratio H:/H,. The flow was fully turbulent. 
L is always defined to the reattachment point, as 
shown in Fig. la. These data were obtained from 
schlieren pictures like those shown on Fig. la. The 
change in the structure of the flow is abrupt and well 
defined, provided the ratio of the boundary-layer 
thickness to the separation-stop height, 5/4, is less 
than about 1; when 6 A; is larger than 1 the cavity 
opens suddenly, but closes gradually. There is a 
small hysteresis region, which is indicated, such that 
if L is being decreased continually, the wake will re- 
main closed longer, and conversely. In the closed- 
open transition the cavity is quite unstable. 

The straight-line correlation of the critical closure 
distances is explained by the simple hypothesis that 
the cavity opens when the vertices of the two separate 
(and initially noninteracting) wakes merge and back- 
flow occurs from the recompression to the separation 
It follows that: 


hip Lp Lr (#) 
- + (2) 
HA, Hi, H, \iN, 


wakes. 


The intercept and the slope of the critical-distance 
line depend only on the geometries of the separation 
and the recompression wakes, respectively. Accord- 
ingly, a review of the information available on the 
geometry of these element-wakes should provide a 
tramework for checking and extending the cavity- 
closure measurements of Fig. 2.T 


(a) The Separation Wake 

The pressure in the wake (base pressure) is related 
to its length/height ratio by a Prandtl-Meyer ex- 
pansion. In terms of the Chapman-Korst model, the 
wake pressure is also related to the pressure after 
recompression through the mixing parameter (ux) 
and the free-stream Mach number. The existing solu- 
tion for the development of the shéar layer?~‘ predicts 
that ux is basically constant in either the totally laminar 
or totally turbulent regime. For constant ws the 
decrease in base pressure and the change in Prandti- 
Meyer angle are related in such a way that the ge- 
ometry of the wake remains approximately constant 
over the midsupersonic Mach-number range of these 
tests. L,/H, is then very nearly constant with both 
M and Ke within either regime, as shown in Fig. 3. 
The experimental points verify this reasonably well 
in the turbulent flow. 

In the laminar region there is considerably more 
scatter among the few experimental results that are 
available.**?7 However, in this regime one has more 
confidence in the analytical model. Laminar-tur- 
bulent transition is defined empirically. One notes 


+ Note that this mechanism is a direct explanation of the 
“opening” of the cavity, but not its closing. When L goes 
through L,, from below, the wake becomes unstable in the sense 
that it jumps over to the closed-cavity configuration readily 
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the well-known difference between two-dimensional Reynolds number.T However, L;/H2 is found to be a 


and axisymmetric laminar flow, for which there is no 
clear analytical explanation. * 


(b) The Recompression Wake 


cannot be applied 


However, there is 


The Chapman-Korst model 
directly to the recompression wake. 
a strong similarity in the characters of the recompres- 
sion and the separation wakes. 

The extensive investigations of Chapman ef a/.' 
show the existence of a “‘plateau’’ pressure Pp over 
the central portion of the roughly triangular wake 
(as shown schematically in Fig. la). If Pp is identified 
with the pressure in the external flow after turning 
by the oblique shock which emanates from the apex 
of the triangle. the length-to-height ratio of the wake 
can be deduced. With the use of Figs. 34 and 42 of 
reference 4, Fig. 4 has been constructed; direct 
measurements of L,///2 from schlieren pictures indicate 
that thisis correct. The range in which transition oc- 
curs is defined approximately by a few specific tests. It 
should be borne in mind that the entire curve is based 
on a dense sequence of experimental points. 

It is evident that Ly/H2 is also approximately 
independent of MM in the midsupersonic range for 
turbulent wake flow and also quite constant with 


* The equivalent flat-plate length Reynolds number used in 
Fig. 3 is the simplest parameter governing the transition region, 
although obviously it is one that does not describe it completely. 


much stronger function of both 7 and Xe in the laminar 
region than is Lp/H,. It is also interesting to note the 
limit of applicability of this model, which is associated 
with the minimum Mach number for the existence of 
an oblique shock, as indicated in Fig. 4. 


The mechanism which establishes the structure of 


the recompression wake involves, apparently, two 
distinct phenomena: 
(a) A pressure-rise-induced separation of the on- 


coming boundary layer (‘free interaction’”’ in Chap- 
man’s terminology) which sets the maximum value of 
the pressure at the apex of the cavity, P 4. 

(b) A further pressure rise from this separation 
‘mean’’ wake pressure Pp, which 


pressure P4 to the 
governs the deflection of the external flow and the 
shape of the wake. 

The “ability to recompress,” 
essential role in the Chapman-Korst model (there 
Pr’'-Pz replaces Pp-P4) seems to be a fundamental 
It cannot be explained, for 
the equilibrium 


which also plays an 


characteristic of wakes. 

recompression wakes, in 
between the dynamic pressure of the separating stream- 
tube in the shear layer and the pressure in the external 
flow, as suggested for separation wakes. However, 
the empirical evidence of Reference 4 indicates that 


terms of 


t The slight slope is due to the VC dependence suggested in 
reference 4. However, this is not certain from the measure- 
ments, which could be considered independent of Re 


within 


experimental scatter. 
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one can write: 


(3) 


(Pp — Pa4)/Po = $(Re,Mo) 


A simple correlation of the data in the range of 
1.2 < Mo < 4 gives ¢ = K(M,—1) where K ~ 0.05 
for laminar flow (Fig. 38, reference 4) and 0.25 for 
turbulent flow (Fig. 46, reference 4). Other cor- 
relations have been suggested in references 20 and 21. 
However, we find their phenomenological basis ques- 
tionable. In our opinion, Eq. (3) is an expression of the 
ability of the internal flow (whatever its nature) to 
support a pressure rise. It has no relation to the free- 
interaction problem which is implicitly involved in the 
form of the correlation proposed in these references. 
We have no model for the above form of ¢ at this time; 
itis a frankly empirical equation. 


(c) Correlation with Cavity Closure 

The data assembled in Figs. 3 and 4 show that at 
least the gross trends exhibited by the proposed notch- 
closure mechanism, Eq. (2), verified: (a) the 
numerical values of both the intercept and the slope of 
Fig. 2 check the values of Lp/H, and Ly/ H2 from Figs. 
4, respectively; (b) the Mach-number inde- 
pendence within the The 
numerical check is actually remarkably good in view of 
the fact that the “lengths’’ of the element-wakes are 
probably not definable rigorously and some inter- 
action before actual opening of the cavity is to be 
expected. The measurements fall all within the 
fully turbulent region, where the wake geometry is not 


are 


> 


3 and 


test range is verified. 
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sensitive to either 1/0 or Re. The study of the “‘free 
wake’”’ discussed below extends these data somewhat; 
however, until more tests are available, Eq. (2), in 
conjunction with Figs. 3 and 4, appears to be adequate 
to serve for estimates of critical distance in the laminar 


region also. 


Free Cavities 

We denote a cavity in which the floor is substituted 
by a plane of symmetry of the wake flow a free cavity. 
Experimentally, such a configuration was obtained by 
introducing a cylinder into the mixing region behind 
the recompression point of a blunt-base wake. Fig. 5 
shows schilieren pictures of the closed and the open 
type of cavity which occurs in this case, as it does in 
the case of the 
critical closure distances. 
that the linear behavior, Eq. (2), found for the bounded 
cavity holds also for the free cavity and that the 
This suggests 


notch. Fig. 6 shows the measured 


It is immediately apparent 


constants are very nearly the same. 
that, if our mechanism of closure is correct, there must 


form in both cases, prior to closure, ‘‘wakes’’ behind 


the base and ahead of the recompression surface 
having (at the point of closure) length-to-height 
ratios identical with those of the corresponding 


“bounded” wakes. 

The correspondence of the intercepts can be ex- 
plained without difficulty because they depend only 
on the lengths of the separation wakes, which are 
approximately equal for the step and the free wake 
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30r — | | For comparison, Fig. 9 shows similar isobars of the 3 
| | Rayleigh pressure (P72) in an undisturbed wake. * 
41 weoce | Ho! is We shall now reason backwards from the hypothesis 
2 } 4 ‘ 3 * > 
SS —— +H , i (aa 97 outlined in the preceding paragraph. Consider the 
yi ; | sé ” —_ . - i “ — “,* 2 . 
ae | den Pa | closed”’ free wake just before the critical point. 
| If the Ly/H, of the free and the bounded recompres- 
= N M,* 278 ; ; ? ; 
va T Pa ccoctioeel” sion wakes are identical, then so are the deflections 
c é x = 
. oo Ler | imposed on the external flow,f the ‘‘external’’ pressure, : 
CR s t 
sre wee L/H 60 | and the mean wake pressure Pp. 
' Reg 1.6 X10° —————— 9\;: ¢ 
meen — If Eq. (8) indeed expresses a fundamental charac- 
L/M,* 69 ec | | teristic of wakes, then the pressures at the apex of the 
M,* 186 ; : r 
Pe “4 . recompression wakes, P?4, are also equal in the two 
10 = Se Re3x10° LyH=45 —— I faci P. j ; i 
oeve L /H,*60 reid n a acing steps )' a 1s set by the 
ne | Res 4X 10° free-interaction”’ mechanism, that is, by the pressure 
oL | | | | rise that the boundary layer will support before 
separating. By what is it set in the free-wake case? i“ 
| | | 
* As measured by a total-head tube. Static-pressure measure- 
° 1 ments were also attempted. They showed a transverse gradient 
ad sl sa ‘in -” 2s 30 due, probably, to disturbances introduced by the probe, a 
4 ; : ae : problem encountered often in such surveys in highly disturbed 
Fic. 6. Critical closure distance for free cavities in the transi- regions. Incidentally, when the Rayleigh pressure is related to 
tional regime. ees ‘ : 
the calculated static pressure, following the wave pattern shown 
in Fig. 9, a remarkably consistent external flow field can be 
constructed. f°) 
+ The free-recompression wake is long enough just before wake 
opening to be approximated by a triangular region, at least | 
re: : reasonably far from the apex. 
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S act 
The principal and interesting difference between the - pre 
free and the bounded cavities lies in the character of a pre 
the recompression wake. lay 
As the cylinder penetrates into the mixing region 
behind the throat of the separation wake from down- re: 
stream infinity, the classical stagnation-point flow ex] 
must change to a vortical structure characteristic of the 
“separations.’’ The bulged-out portion of the bow the 
shock visible in Fig. 5 does contain a pair of vortices ( 
this can be deduced from the fact that the centerline wa 
pressure on the surface of the cylinder is less than the to 
axi surfe 255 aw: 7) > centerline, : ? 
maximum surface pressure away from the wae rline = 68 CS CM “30 Ps) 
as shown by typical curves in Fig. 7. Fig. 8 shows a ° the 
complete isobar map of the surface pressure on the =~ _ ee ( 
x . Me Fic. 8. Isobar map of the surface pressure on the recompression 
cylinder on both sides of the critical closure length. cylinder, for varying L/H. My = 2.8, Re, = 0.5 X 108. of 1 
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Isobar map of the Rayleigh"pressure ratio in the free- 


Fic. 9. 
(The ratio is referred to the 


mixing region and the shear layer. 
local free-stream pressure. ) 


The velocity at the centerline of the free-mixing 
The flow can support a streamwise 
(which is assumed to 


region is not zero. 
pressure rise to stagnation 
dominate over the pressure rise due to transverse 
indeed, the flow must stop 
Furthermore, 


inomentum transfer); 
before reversing to open the cavity. 
there can be no backflow before the apex of the re- 
compression wake penetrates into the subsonic flow in 
the mixing region, or at least to the sonic “throat.” 
It is plausible to assume that the cavity opens when 
this happens. Now, at the “sonic’”’ point the center- 
line pressure rise to stagnation is: 


(P, — Pr')/Pr’' » (y/2)Me?~0.6 at Men 


This phenomenon takes the place of the “‘free inter- 
action.’’ By coincidence, the numerical value of this 
pressure rise is approximately equal to the value of the 
pressure rise to separation of a turbulent boundary 
layer, although the mechanisms are quite different. 

The foregoing arguments must be taken 
reserve; nonetheless, they do provide a_ plausible 
explanation for the surprising similarity in results for 
the closure of two quite different cavities, and moreover 


with 


they suggest that: 

(a) The observed equality of the free- and bounded- 
wake critical distances is coincidental and restricted 
to this particular Reynolds-number range. (P4 — 
Ps) is fixed for free wakes while it depends on +/C;, for 
the upstream-facing step. 

(b) Eq. (3) may indeed be a fundamental description 


of recompression wakes. 


FLOWS—PRESSURE FIELD 463 


The tree-wake tests also verified the proposed 
mechanism of closure in these aspects: (a) the critical 
distance is independent of the Mach number within the 
test range; (b) by changing the chord of the wedge it 
was possible to vary the Reynolds number at separation 
through the separation-wake transition region, and it 
was found that the intercepts are as required by Eq. 
(2) (they are cross-plotted on Fig. 3). In drawing 
the correlations (Fig. 6), their slope was kept constant 
on the ground that the highly disturbed flow in the 
mixing region would certainly be turbulent in spite of 
the fact that the separation wake may be laminar. 
In this sense the Reynolds number at separation is 
not an adequate criterion to describe the transition 
region of the cavity-closure phenomenon. More data 
would be needed to clarify this point. 


Axisymmetric Cavities 


Fig. 2 includes two points, obtained during this 
work, relating to the closure of axisymmetric cavities 
formed between a conical spike and a blunt afterbody. 
When correlated in terms of the pseudo-two-dimen- 
sional step height between the outer radius of the body 
and the sting, the axisymmetric cavity closes at the 
same critical separation as the two-dimensional notch. 
Other data’: both for notches in cylindrical surfaces 
and for spike extensions verify this observation for an 


extended range of sting-to-body radii. However, when 
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extended to very large ratios of separation- to recom- 
pression-step height (//2///,;), axisymmetric closure 
exhibits a very pronounced hysteresis. For instance, 
Beastell et a/."° recorded a hysteresis range of AL //1, = 
12 at H2/H, = 20. There are no equivalent data for 
two-dimensional steps. 

A limited amount of pressure measurements in 
axisymmetric cavities was also collected. They are 
not extensive enough to allow a conclusive comparison 
with notches to be made, and will not be shown in 
what follows. In general, the trends observed in two- 
dimensional and in axisymmetric flows seem to be the 
same. However, the rise of pressure along the floor of 
the open axisymmetric cavity (the sting) from the 
separation to the recompression step seems to be 
considerably larger, as evidenced also by a comparison 
of Johannsen’s data!’ with, for instance, Morozov’s 
data’ at approximately the same Mach number. 


The Distribution of Pressure in the Cavity 


The many pressure-distribution records obtained 
fall into a few distinct types of which we shall show only 
examples. This investigation was restricted to notches 
of equal separation- and recompression-step height 
(H; = H;) and to a fully turbulent oncoming Loundary 
layer. There are four principal parameters the effect 
of which is to be studied, namely: the length-to- 
height ratio of the notch, the free-stream Mach number, 
the geometry of the recompression face, and the ratio 
of boundary-layer thickness to notch height. The 
last parameter replaces the Reynolds number which is 
not as significant within the turbulent flow region. 

The first group of figures shows simultaneously the 
typical results and the effects of recompression-face 
geometry and boundary-layer thickness. These will 
be discussed in sequence. The effect of Mach number 
will be shown in another group of curves. 


(a) Floor Pressure 

Typical floor-pressure distributions which charac- 
terize a small, intermediate, and near-critical notch 
are shown in Fig. 10, for two values of 6/H and the 
three geometries tested (to be discussed subsequently 
in detail). These trends recur for all Mach numbers 
even if the magnitudes of the pressure ratios differ. 
The pressure in the notch is quite uniform for small 
notches, with a slight minimum near the center of the 
span. As L increases, the minimum moves toward the 
separation step, and the pressure gradient to recom- 
pression strengthens. Near L,, this gradient becomes 
very strong and a “plateau’’ in pressure, which is 
characteristic of the recompression wake, is often 
noticed. These trends are shown in more detail 
by the “isobar map” of Fig. 11. The variation of the 
floor-pressure distribution is gradual and continuous 
throughout the region of open wakes. 

The transition from the ‘‘small’’ to the intermediate 
length cavity is associated with the appearance of a 
weak oblique shock rooted in the shear layer, which is 
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ned Fic. 18. The thickness of the recompression profile as function of L/H, Mach number, and 6/H. 
nly 
hes 
ght visible in Figs. la and 14. A similar shock was’ noted 
ary also by other investigators.'*'* This shock is not 
ect caused by or dependent on laminar-turbulent transition 
‘to- in the shear layer although a similar phenomenon may 
er, be due to this;* instead it is due to a deflection of the 
tio external stream, which is a part of the recompression 
rhe mechanism. 
1 is 
(b) Recompression Step Pressure 
the Fig. 12 shows typical distributions of pressure along 
— the recompression face (normal to the flow). The 
vill pressure profile is fairly uniform near the floor but 
er exhibits a gradient at the outer edge of the step which 
suggests a recovery of dynamic pressure of the portion 
of the free-stream layer which impinges on the step 
(and which is returned into the cavity). Particularly 
”" striking is the decrease of the ratio of maximum to 
™ minimum pressure with increasing notch span. This 
- will be discussed later. The thickness of this profile 
ly (A) can be clearly defined, as shown graphically in 
“ Fig. 12. Fig. 13 depicts the variation of A as a function 
©, of L/H (that is, the length of the free-shear layer), for 
ill several values of the external flow parameters. 
he 
he (c) Geometric Similarity Parameter 
a- The result of the pressure distribution studies shows 
eS conclusively that the appropriate similarity parameter 
is governing the pressure field at a fixed Mach number is 
n the ratio of the length of the notch to its critical 
il length, L/L,,. This holds also for notches with 
le external-flow compression or expansion at separation 
IS (1, # H2). When 6/H is large one would expect 
secondary effects which do not scale with L/L,,. 
e However, within the range of 6/H investigated, sim- 
‘ ae ge Peto va. (eke et when Fic. 14, The ee Wd — boundary layer 
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Fic. 15. Effect of 6/H on the pressure distribution in a separation wake (M = 2.8, Re, = 1.5 X 105) 


(d) Effect of the Thickness of the Boundary Layer 


Fig. 14 shows a series of photographs of the closure 
of a cavity in a very thick boundary layer (6/H ~ 2). 
There is no fundamental difference in the flow pattern 
but the closure is less abrupt. In general, an increase 
in the 6/H ratio brings about the smoothing out of all 
pressure gradients. It can be presumed that this 
follows from the decreased momentum transfer to the 
i.e., a decrease in the value of the mixing param- 





cavity 
eter Ux. 

For example, the pressure variation on the floor 
behind the downstream-facing step is shown in Fig. 15 
for two values of 6/H. The decrease of all pressure 
gradients, and the increase in base pressure with 
increasing 6/H, are evident. Typical measurements of 
the base pressure for various 6/H are indicated in an 
insert.— They show the well-known increase in 
Pz/PpR with 6/H, but also that it is not strong until 
5/H becomes larger than unity. 

Similar trends were observed for notches. 
pressure distributions were already shown in Fig. 10. 
Note that for 6/H >1 the pressure in the cavity 
increased above the free-stream static pressure for 
small notches (see also Fig. 16). This rather surprising 
result was also recorded (but not stated) by Thomann” 
(MU = 1.8, 6/H = 1 to 2, turbulent flow). It seems 
that the subsonic portion of the thick boundary layer 
diffuses the high pressure near recompression beyond 
the boundaries of the notch proper and leads to a 
slight ‘‘bulge’’ in the supersonic external flow. 

Fig. 12 shows that the effect of 6/H on the recom- 
pression-face pressure distribution is quite small, 


Typical 


{ Tests of both constant H, varying 6 and constant 6, varying 
H were conducted, and they proved that the ratio 6/H is the 
proper similarity parameter. 


even over a large range of this ratio. The explanation 
for this lies in that the recompression flow is highly un- 
steady and that these measurements are time averages 


which efface the detail of the shear-layer profile. 


(e) Effect of Notch Geometry 

The evidence indicates that the geometry of the 
recompression step does not essentially modify the 
pressure distribution in the cavity (Fig. 10). This is in 
accord with the results of several other investigations 
(e.g., reference 4). This shows that the length Z of 
the cavity should be defined as shown in Fig. 1b—1.e., 
to the reattachment of the free-streamline. 


(f) Effect of Mach Number 


The qualitative trends of the pressure distributions 
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Fic. 16. Variation of the pressure distribution in the notch with 


Mach number (Type 1, 1.6< 6/H< 2.3). 








be 


lai 


efl 
th 


da 


sti 
by 


ra’ 
wi 


flo 
Or 
th 
fac 
le 
th 
th: 
sti 
po 
ris 
un 


the 
in 
aft 
shi 
pa 
flo 
mi 
rec 
no 
the 
bo 
Th 
of 
ou 
the 
en 


lor 


wa 


ion 
an- 
Fes 


che 
he 

in 
ns 


of 
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repeat themselves for each Mach number. The 
relative magnitudes of the pressure ratios are not, how- 
ever, simple functions of it. 

To demonstrate the Mach-number trends clearly, 
the variation of the pressure just behind separation 
(Y L = 0) and just ahead of recompression (X /L = 1) 
is shown in Fig. 16 as a function of Mach number and 
notch length /height ratio. 6/H was held as constant 
as possible. These tests were all conducted in a thick 
boundary layer; note that the floor-pressure ratios are 
larger than unity for small L/H, which is characteristic 
of the flow with 6/H >1, as indicated above. This 
effect seems to be stronger as Mp» increases, although 
this cannot be concluded with certainty from this 
data because, unfortunately, 6/H increases then also. 
The wakes close at L/H~12 independently of Mo, as 
stated earlier, which is evidenced in this representation 
by the drop of the separation-pressure line. Of 
course, the final (closed cavity) separation-pressure 
ratio (base pressure) is lower for higher 1%, in accord 
with theory. 


Structure of the Flow Over a Notch 


Certain conclusions concerning the structure of the 
flow over a notch emerge from these measurements. 
One observes that the pressure rise along the floor of 
the notch and the pressure rise along the recompression 
face exhibit opposing trends with increasing notch 
length. This is shown in Fig. 17. In short notches, 
the floor pressure is nearly uniform (sometimes larger 
than free-stream when 6/H is large) while there is a 
strong impact-pressure recovery near the reattachment 
point. In long notches there is a strong floor-pressure 
rise, but the recompression-face pressure is nearly 
uniform (see also Fig. 12). 

This suggests that there are two possible structures of 
the flow over an open notch cavity, which the sketches 
in Fig. 17 describe better than words. In long notches, 
after an initial expansion into the cavity, the free- 
along a curving 


” 


shear layer gradually ‘‘recompresses 
path. The increasing pressure is reflected along the 
floor and the momentum of the inner portion of the 
mixing layer is absorbed, resulting in no impact 
recovery against the recompression wall. In short 
notches, if there 1s any boundary layer before separation, 
the impact-pressure rise at reattachment thickens the 
boundary layer and diffuses the pressure gradients. 
The cavity pressure is nearly uniform. The thickening 
of the boundary layer can be sufficient to cause an 
outward predeflection of the external flow and raise 
the pressure level above free-stream pressure in the 
entire cavity. 

More detail concerning the structure of the flow in 
long notches is discussed in the next paragraph. 


Pressure Traverses Through the Shear Layer 


The free-shear layer in a large notch (L/L,, = 0.9) 
was investigated (at 1/4) = 2.9) in more detail. A 
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series of stagnation-pressure traverses normal to the 
flow, beginning with the boundary layer ahead of 
separation, and a streamwise traverse of both the 
total and the static pressure (parallel to the notch 
along plane z) are shown in Fig. 18. 

It was found from the streamwise traverse that the 
static pressure gradually decreases at first, as the 
flow passes through a nearly centered expansion fan 
emanating from the separation point, and then in- 
creases through the central predeflection shock. This 
is consistent with schlieren observations and with the 
measurements of floor pressure. From these measure- 
ments, the free streamline shown in Fig. 18 was con- 
structed. * 

The transverse surveys at stations A through E are 
also presented in Fig. 18 in the form of local Mach 
numbers. The average static pressure at plane z 
and on the notch floor (full circles) as well as the 
floor pressure alone (flagged circles) were used to 
reduce the data. 

Next, integrating the mass flow piecewise at each 
* The measured static pressure rise between stations C and D 
was gradual, due to the instability of the shock position. In 
drawing the streamline this gradual rise was disregarded and 
instead the full static-pressure rise between these stations was 


used to compute the flow-deflection angle 


SCHEMATIC STRUCTURE OF THE FLOW: 
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OVERALL PRESSURE RISE IN NOTCH: 
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Fic. 17. The recompression process as function of notch L/H 
and Mach number (Type 1 notch, 6/H = 1.6 to 2.3). 
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Fic. 18. Flow and Mach number distributions through the shear layer and in the cavity (7 = 2.8, 6/H = 0.4, L/H = 10, Type 2) 


station from the free streamline, starting with station 
A (assuming isoenergetic compressible flow), the 
streamtubes in the shear layer were found, as shown 
in the figure. With due consideration for the in- 
accuracies involved in this procedure, which is very 
sensitive to the profile in the low velocity region of the 
layer, some interesting observations can be pointed 
out: 

(a) Overall mass conservation between stations A 
and E is satisfied well within the accuracy of the integra- 
tion procedure. 

(b) However, at all internal stations (B, C, and D) 
there is an unaccounted-for mass flux in the direction of 
the external flow. A downstream-pointing total tube 
was then used to verify the existence of a reverse-flow 
layer. Such a reverse flow was indeed found to exist at 
the central stations C and D, although the mass 
flux it carries does not compensate for the streamwise 
flux at these stations. At station B no reverse velocity 
was measured within the accuracy of the readings. 

(c) The inner streamtube is found to bend into the 
notch very strongly; one cannot rely on these com- 
putatious to give its shape but the fact that it does 
penetrate into the cavity seems real. 

(d) At station E, which is after recompression, the 
profile shows a definite and large velocity all the way to 
the wallf and does not resemble that of an established 
boundary layer. 

The results of these pressure surveys are satisfying 
as regards the free stream above the notch and the 
flow in the outer portion of the shear layer and the 
consistency with floor pressure surveys. However, they 
do not produce a satisfactory picture of the inner flow. 
It will be shown in Part II of this series that the 
motion in the cavity is complex and unsteady so that 
mean-pressure measurements, which are direction- 


{t Measurements were made with a flattened boundary-layer 
probe with an 0.012-in. outside dimension. 


and pulse-shape sensitive, cannot alone describe it 
adequately. 


The Drag of Notched Walls 


The drag of a notched plate can be calculated from 
the measurements of pressure on the faces normal to the 
flow (the shear on the floor of the notch is negligible). 

It is convenient to use a coefficient 


_ 218 ft Ppe-P,.fY | 
Ce wt : i} d (4) 
“i 4 M,? zZ 0 Po H 


The floor pressure at the separation corner was 
assumed to be the pressure over the separation step, 
as shown by several investigations. The integral over 
the recompression step can be obtained from data at 
hand (such as Fig. 12). 
necessary because the subtraction of two nearly equal 
numbers (Pr/Po)av. and (Ps/Po) for small notches 
Fig. 19 
shows, then, the drag coefficient Cp as a function of 
L/H and Mo, based on the smoothed-out curves. The 
data are for a Type I notch but, as indicated previously, 


Some smoothing out was 


emphasized unduly the experimental errors. 


the geometry of the recompression face does not 
introduce any essential changes. 

It is interesting to note that the drag coefficient 
appears to decrease for long notches but increase for 
small notches with increasing Mach number. Consider 
this in terms of structure of the external flow over the 
notch outlined in a previous paragraph. The flow over 
a long notch is characterized by inward deflection of the 
free streamline as shown schematically in Fig. 17. 
The potential pressure exerted on this double-wedge 
inflection in the free boundary (and transmitted to the 
cavity) is proportional to Mp o?(Mo?—1)~—? 8%. Ob- 
servations show that the deflection angle 8 depends 
principally on L/L,,. Thus, the pressure drag coef- 
ficient (computed along the free streamline) would be 
of the order of 
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This is in qualitative agreement with the variation of 
long-notch drag based on pressure readings along the 
boundaries of the notch. 

In small notches, on the other hand, particularly 
when 6/H and therefore 6/L are large, which is the 
case here, the pressure in the cavity is governed by the 
outward deflection of the free stream due to the thicken- 
ing of the viscous layer. This thickening can be as- 
sumed to increase with Mach number, resulting in an 
increase of the flow-deflection angle 8(~A6/L) with 
My. Uf this functional dependence is stronger than the 
one-half power, which seems reasonable, a positive 
slope of Cp versus / will result. Note that this 
effect must be related to the initial magnitude of 
6/H, and when this ratio is small the observed increase 
of Cp with A/ in small notches may not occur. 


Subsonic Flow over a Notch 


In a subsonic flow, there is no abrupt change from 
closed to open configuration of the cavity (see also 
reference 17). However, the changes in the floor- 
pressure distribution are very similar to those observed 
in supersonic flow, as seen from Fig. 20. One dis- 
tinguishes, although less clearly, the characteristic of 
the small-cavity structure, in which the pressure is 
nearly uniform within the entire notch. Longer 
notches (L/H >4) exhibit a pronounced pressure 
distribution which, in this subsonic case, goes over 
smoothly to the distribution characteristic of a ‘“‘closed”’ 
cavity—i.e., a flo wwhere one can distinguish clearly the 
two independent separated regions. 

The base pressure at the separation tends then to 
a negative value which checks the prediction based on 
Korst’s analysis’ quite closely. The basic similarity 
criteria that were found to generalize the supersonic 
data appear to apply also to the subsonic flow except 
that, as already noted, the boundaries of the various 


cavity regimes are less clearly defined. 


Concluding Remarks 


The results of this investigation establish conditions 
for the occurrence of separated cavities in depressions 
in a boundary and define the parameters which govern 
their structure and the pressure field. The principal 
similarity parameter is found to be the ratio the 
length to the critical closure length of the cavity. The 
latter is related to the properties of simple wakes which, 
in turn, can be described by semiempirical models and 
correlations of data. However, detailed surveys of 
the shear layer show that models which postulate a 
steady mass convervation in the cavity and seek the 
conditions governing the equilibrium of the flow 
in the profile of the shear layer alone are not adequate 
to describe cavities in which the flow recompresses 


against a step. 
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Similarities and differences among various types of 
cavities and wakes are studied and discussed. The 
effect of the initial boundary layer is investigated and 
it is shown that it begins to influence the structure of 
the flow strongly only when its thickness exceeds the 
depth of the notch. 

This investigation was continued and further results 
will be described in Part II of this series. These concern 
the flow inside the cavity and explain the results of the 
pressure traverses reported in this text and confirm 
the pulsating nature of cavity flow. Measurements of 
heat transfer across the separated region will also be 
reported. 

From the comparison of such diversified but comple- 
mentary results, one hopes to gain insight into the 
complex structure of separated flows. 
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Fic. 19. The drag coefficient of a notched wall (Type 1) 
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Fic. 20. A map of the floor pressure distribution in a rectangular 
notch in subsonic flow. 
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The Turbulent Boundary Layer on 
Chemically Active Ablating Surfaces’ 


M. RICHARD DENISON* 
Aeronutronic, A Division of Ford Motor Company 


Summary 


Incompressible turbulent-boundary-layer analysis is extrap- 
olated analytically to the case of a compressible turbulent 
boundary layer with ablation or mass injection at the surface. 
The effects of chemical reactions such as dissociation and re- 
combination as well as combustion are included. The analysis 
applies to blunt as well as sharp bodies which are either axi- 
symmetric or two-dimensional. When the turbulent Lewis and 
Prandtl numbers are unity, it is found that, as in the laminar 
case, little detailed knowledge of the chemistry inside the bound- 
ary layer is required in most instances. The conditions at the 
surface and the outer edge of the boundary layer are often suffi- 
cient for prediction of heat and mass transfer 

Comparison is made with experiments on the combustion of 
graphite under turbulent flow conditions. Prediction of ablation 
rates within about 30 percent accuracy is obtained when em- 
pirical constants obtained from incompressible velocity profiles 


with no mass injection are used. 


Symbols 


a = parameter in enthalpy-velocity relationship, see Eq. 
(32) 

Qy = skin-friction parameter, see Eq. (16) 

b = parameter in enthalpy-velocity relationship, see Eq. 
(33) 

B = ratio of mass injection to skin friction, see Eq. (8) 

B* = value of B with reactants in stoichiometric proportion, 
see Eq. (9) 

B = ratio of mass injection to zero-blowing skin friction, 
see Eq. (40) 

Cye = Skin-friction coefficient based on local free-stream 
conditions 

iis = specific heat 

D = molecular-diffusion coefficient 

Dy; = turbulent-diffusion coefficient, see Eq. (4) 

h = enthalpy of mixture 

h., = enthalpy of condensed phase at surface 

h; = enthalpy of chemical species i including heat of 
formation 

h® = heat of formation of ith chemical species 

h® = heat of formation of mixture, see Eq. (24) 

h, = total enthalpy of mixture 

H = enthalpy of mixture with zero datum, see Eq. (29) 

k = unity for axisymmetric and zero for two-dimensional 
bodies 

K = thermal conductivity 

K = constant in mixing-length law 

K; = mass fraction of ith molecular species 
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Kr = turbulent thermal conductivity, see Eq. (4) 

l = mixing length 

Ler = turbulent Lewis number, see Eq. (4) 

Mi» = net mass flux of ith species at the surface 

My» = net mass flux of mixture at surface 

M = molecular weight 

a = pressure 

Pry = turbulent Prandtl number, see Eq. (4) 

Ga = heat absorbed by conduction into condensed phase 
and/or radiation away from the surface 

r = ratio of skin-friction coefficients with and without 


mass injection 
’mi = Yatio of mass of atomic species m to unit mass of 
molecular species 7 
Ym¢ = ratio of mass of atomic species m to unit mass of 
condensed-phase material 
distance from axis of symmetry to body surface 


Yo = 

Rewr = effective Reynolds number, see Eq. (21) 

T = temperature 

u = component of velocity in x direction 

U = universal fit for (p/P)! ? vs. enthalpy, see Eq. (24) 

v = component of velocity in y direction 

x = tangential distance from stagnation point along body 
surface 

y = distance normal to body surface 

Z = velocity ratio u/uUe 

Ze = roots of Eq. (34) 

a = coefficient in Eq. (40) 

B = unity for the Prandtl mixing-length assumption 
1 = Ky and zero for the von Karman similarity law 

In, = atomic mass fraction, see Eq. (2) 

I.. = atomic mass fraction of oxygen at outer edge of 
boundary layer 

6 = boundary-layer thickness 

€ = eddy viscosity, see Eq. (4) 

be = molecular viscosity 

v = kinematic viscosity 

é = parameter for determining atomic mass fractions, see 
Eq. (22) 

p density of mixture 

T = shear stress 

? = integral defined after Eq. (21) 

v = integral defined after Eq. (16) 

Subscripts 

e = conditions at edge of boundary layer 

t = index for ith molecular species 

m = index for mth atomic species 

0 = conditions at no mass injection 

w = surface conditions 

Superscripts 
* = stoichiometric conditions 
Introduction 


4 i NHE STATE OF DEVELOPMENT of turbulent-boundary- 
layer analysis is often deplored by fluid dynami- 
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cists. In spite of some advances which have been 
made in understanding turbulent shear flow,'~* much 
remains to be accomplished before this knowledge can 
be incorporated into analysis. Unfortunately, the 
engineer is faced with turbulent flow far more often 
than laminar flow. In the case under consideration, 
involving ablation of chemically active surfaces, this is 
even more likely to be true. The engineer usually 
overcomes the analytical difficulty because turbulent 
skin friction, mass transfer and heat transfer can be 
predicted with some degree of success with the aid of a 
few empirical constants obtained from experiments, 
even though the details of the phenomena remain 
obscure. 

The analysis presented herein is substantially the 
same as that previously given in reference 4. Since 
that time it has been subjected to experimental correla- 
tion in the case of graphite combustion,® and found to 
be capable of predicting ablation rates within 30 per- 
cent accuracy when empirical constants obtained from 
incompressible velocity profiles are used. In view of 
the complexity of the analysis such agreement is con- 
sidered very good and adds to one’s confidence in the 
turbulent analysis. The experiments are described 
briefly in this paper. 

The analysis itself is similar in many respects to the 
theories that have been developed for laminar flow.°~* 
The principal difference is that macroscopic rather than 
molecular transport processes are involved. In the 
gas phase the assumption is made that diffusion and 
convection are rate-controlling, so that any chemical 
reactions that occur can be considered to be in chemical 
equilibrium. At the interface between the gas and the 
condensed phases it is possible to consider nonequilib- 
rium heterogeneous reactions, but in the interest of 
simplicity these reactions are also assumed to be in 
equilibrium. For relatively moderate ablation rates 
and high pressures this assumption is probably ade- 
quate. When the turbulent Lewis and Prandtl num- 
bers are assumed to be unity, it is found that, as in the 
laminar case, little detailed knowledge of the chemistry 
inside the boundary layer is required in most instances. 
Only the conditions at the surface and at the outer 
edge of the boundary layer are important. However, 
in contrast to the laminar case, when gas-phase com- 
bustion reactions occur, the skin-friction law is some- 
what affected 

The work presented represents an analytic extrap- 
olation of semiempirical analyses for incompressible 
turbulent flow to the case of a compressible turbulent 
boundary layer with ablation or mass injection at the 
surface, including the effects of chemical reactions such 
as dissociation and recombination as well as combustion 
in the boundary layer. The analysis applies to blunt 
as well as sharp bodies which are either axisymmetric 
or two-dimensional. The skin-friction law can be 
shown to be approximately equivalent to the solution 
of Dorrance and Dore® for compressible flow with mass 
injection, to reduce to the solution of Van Driest'® 
for no blowing, and finally to the incompressible von 


Kaérman-Schoenherr law when proper thermal and 
geometric properties are inserted. Although this is of 
some comfort, more extensive experimental correlation 
is desirable. 


Fundamental Equations 


Although prediction of turbulent boundary-layer 
phenomena, within the present state of knowledge, 
must be semiempirical in nature, analysis of first prin- 
ciples helps form a logical framework for physical rea- 
soning. The development of the fundamental equa- 
tions rests on the conventional assumption that the 
conservation equations for continuum hold instan- 
taneously. The instantaneous value of each of the 
unknowns is separated into the sum of mean and fluctu- 
ating parts. Introduction of such expressions into the 
conservation equations for reacting gas mixtures and 
averaging according to the Reynolds rules for time 
averages yields average conservation equations that 
include correlations between fluctuating components. 
In the development of these equations it is assumed 
that diffusion due to temperature and pressure gra- 
dients is negligible compared to that due to concentra- 
tion gradients, and that all diffusion coefficients are 
equal. When correlations involving fluctuations in 
molecular-transport terms are neglected and the usual 
boundary-layer assumptions are made, the following 
boundary-layer equations result in the mean steady 
state for axisymmetric and two-dimensional coordi- 
nate systems: 


Conservation of Mass: 
(0/dx) (pure) + (0/dy) (parc) = 0 (1) 


where & is unity for bodies of revolution and zero for 
two-dimensional bodies, and the bars refer to time- 
averaged quantities. 


Conservation of Atomic Species: 


or. or. Oo Ler 4 or, ) 
Uu - pv = (a ) (Z 
ov | Oy) Oy LPrr ® | Oy , 


where I’, is the atomic mass fraction of atomic species 


mM: 
igs = > Vm K; 
i 


with A; the mass fraction of the 7th molecular species 
and r,,; the ratio of the mass of atomic species m to unit 
mass of molecular species 7. Eq. (2) can be obtained 
by forming linear sums of the set of molecular-species 
continuity equations. There are usually fewer equa- 
tions for [’,, than there are variables K;. The remain- 
ing relationships are supplied by the equilibrium con- 
stants for each chemical reaction when equilibrium 
conditions are approximated 
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CHEMICALLY ACTIVE 


Conservation of Energy: 
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where h, is the total enthalpy defined by: 


a7 2 
h =X K, (J Cy dT + as) += 
F 0 


Conservation of Momentum: 


O(t/u,) O(u/Ue) 


a7 pu - 
. Ox oy 


aT. rr O(a, “2 oP 1 (: mt) t) 
(i ) _ _ (2 
oy “ : oy Ox Uu, pple” 
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where Pre = 
Ler = PCy (D a Dr), (K a Ky) 
and 


—(pv)'u’ = € OnW/Oy 


—(pv)' 20K hi = K,oT/dy 


—(pv)'T,,’ = pD,or,,/dv 
where the bars denote averages and the primes denote 
fluctuating components. 

With the assumption that the turbulent Prandtl 
number Pry and Lewis number Le; are unity, and 
neglecting the influence of a pressure gradient on the 
momentum equation, it can be seen that similarity 
exists between the energy equation, momentum equa- 
tion, and the equation for conservation of atomic 
species. Under these conditions it is reasonable that 
h, and IT, be linear functions of the velocity ratio u/1,. 
It is valid to neglect a pressure gradient if the gradient 
is small or if the body surface is strongly cooled, so that 
the coefficient of the pressure gradient in the momen- 
tum equation is small. The advantage of similarity 
is that once the velocity distribution is obtained, the 
use of the von Karman integral equation alone is 
sufficient to estimate the skin friction, heat transfer, 
and mass transfer at the surface. Otherwise, corre- 
sponding integral equations for energy and species 
continuity must also be solved. Another advantage 
is that the heat- and mass-transfer rates become in- 
sensitive to the detailed structure of the boundary 


layer. 
Chemical Reactions, Heat Transfer and Mass 


Transfer 


It can be shown in the general case that the number 
of possible chemical reactions is exactly that required 
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in addition to the atomic mass fractions to completely 
Although 
the number and type of chemical reactions is not 
limited, it is helpful to restrict attention to a specific 


specify all the molecular mass fractions. 


chemical system. The chemical reactions considered 
are a single-combustion reaction, a vaporization proc- 
ess, and the dissociation and recombination of oxygen 
The experiments previously referred 
With graphite 


and nitrogen. 
to involve the combustion of graphite. 
as the surface material, the governing reactions are 


1/3 C3(g) + 0 (g) = CO * 
C (s) = 1/3 Cs (g) 

Oz (g) = 2 O (g) 
N2 (g) = 2 N (g) 


Other molecular forms of gaseous carbon could be in- 
cluded, but recent evidence!! indicates that C; is the 
most prevalent vaporization species at temperatures 
approaching the sublimation temperature. 

In order to determine the conditions at the surface, 
the mass fractions of the various species are required. 
If equilibrium is assumed to exist at the surface as 
well as in the gas phase, the mass fractions of the 
molecular species C;, CO, O, Ov, N, and N» at the sur- 
face can be related to the surface temperature by means 
of equilibrium constants for the reactions indicated 
in Eq. (5). However, these represent only four equa- 
tions for the six unknown mass fractions and surface 
temperature. Further relationships are obtained from 
a consideration of surface boundary conditions. 

The boundary condition for mass transfer at the sur- 
face can be stated in terms of conservation of atomic 


species: 


ring Mry 2 Ymil inv (6) 
t 

where m, is the net flux of surface material, 7,,; is the 
mass of atomic species m per unit mass of surface mate- 
rial, and m;, is the net flux of molecular species 7 at the 
surface. If the Lewis number is unity and the atomic 
mass fractions are linear functions of the velocity ratio, 
then Eqs. (6) lead to the following set of expressions for 
atomic mass fractions at the surface: 


Pine = (U'me + ?msB)/(1 + B) (7) 


> 


where B = m,,/(pettCs,/2) (8) 


A further relationship between skin-friction coefficient 
cy, and B is obtained from the momentum equation, 
which will be discussed later. The atomic mass frac- 
tions at the outer edge of the boundary layer, T,,,,, are 
In the present case there are three atomic 
Therefore, with the defi- 


known. 
species, viz., C, O, and N. 
nition of the atomic mass fractions at the surface, 
I'nw, Eq. (7) represents three relations among the mo- 
lecular mass fractions A, and the parameter B. In- 
cluding the chemical-equilibrium relations this makes 
a total of seven equations for the six molecular species, 
surface temperature, and B. If any one of these, say, 
the surface temperature, is given, all the other quanti- 
ties can be determined. This result can be generalized 
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for any chemical system, and nonequilibrium surface 
kinetics could also be included. 

It so happens that for a wide range of surface tem- 
peratures of interest, the molecular mass fractions of 
the reactants, carbon and oxygen, are negligibly small 
at the surface, corresponding to a condition of complete 
surface combustion. Under these conditions it can 
be determined from Eq. (7) that 


B = B* = (3/4)T,, (9) 


where B* is the value of B with reactants in stoichi- 
ometric proportion. At extremely high surface tem- 
peratures, in excess of 7,000°F or so, the value of B in- 
creases above that given by Eq. (9), corresponding to 
the condition of gas-phase combustion. The tempera- 
ture at which this occurs is somewhat dependent upon 
pressure and T"p,. 

In order to complete the description of the surface 
conditions, the surface temperature is specified by 
means of heat-transfer considerations. The heat 
transfer depends upon the skin-friction coefficient 
through Reynolds’ analogy. Once the skin-friction 
coefficient is determined, the mass-transfer or ablation 
rate is also determined from Eq. (8). 

The heat conducted into or radiated away from the 
solid surface is obtained by a heat balance at the sur- 


face: 


or 
Ga = (x =a - > M ih iw + Mules (10) 


The heat flux is due to conduction out of the boundary 
layer less the difference between the heat convected 
by all individual species and that convected from the 
condensed phase whose surface enthalpy is /,,. If 
Pr, = Ley = 1, and the total enthalpy h, is a linear 
function of the velocity ratio, it can be shown that 


qa = Pell (Cy 2) [Ase cand Nsw i B(hsw _ h-s)| (11) 


When there is no mass transfer from the surface, B 
becomes zero, so that the heat-transfer expression re- 
duces to that of compressible flow with a Prandtl 
number of unity. Eq. (11) can also be put in the 
following useful form: 


Ga = Pelle (C Fe, 2) [hse + An 06Q pow -— BX» = New| (12) 


where Qo. is the heat of combustion for gaseous reac- 
tants at surface conditions, a, is the fraction of heat 
of combustion actually produced (usually equal to 
unity when combustion occurs), A,» is the heat of 
vaporization at surface conditions, and h,, is the en- 
thalpy of gas which exists at the outer edge of the 
boundary layer evaluated at surface conditions. When 
the heat flux given by Eq. (11) or Eq. (12) is matched 
to expressions for radiation at the surface and transient 
heat transfer in the condensed phase, the surface 
temperature can be determined. Because the enthalpy 
hs» depends upon the mass fractions at the surface as 
well as on temperature, it is necessary, except in special 
cases, to solve Eq. (11) simultaneously with Eq. 
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(7) and the equilibrium relations for reactions like those 


of Eq. (5). 

The equations developed above apply to both 
laminar and turbulent flow. The difference between 
the two flow regimes appears in the skin-friction co 
efficient. The result for laminar flow is presented in 
reference 6. An expression for the turbulent skin 
friction is derived in the next section. 


Estimation of Turbulent Skin Friction 


In laminar flow the momentum equation can be 
solved exactly.6~-* For turbulent flow the von Karman 
momentum integral is utilized to determine the skin 
friction. The momentum integral, modified to include 
the effect of mass injection, can be obtained either by 
integration of the momentum equation across the 
boundary layer, neglecting double and triple correla- 
tions, or by consideration of momentum conservation 
for an element of fluid. The result is 


d Fi ro 6 
es puro’ (1 — u/u,)dy = 
dx Jo 


dP f pu’ 

dx J0 (= 
It is consistent with previous assumptions to neglect 
the first term on the right-hand side of Eq. (13). 

In order to solve the momentum-integral equation 
it is necessary to assume a reasonable velocity dis- 
tribution. An approach is utilized which is similar to 
that of Van Driest,'’ who allowed for density variations 
in the compressible boundary layer in such a manner 
that the results reduce to the von Karman—Schoenherr 
relationship for incompressible skin friction with no 
blowing. The increase in shear near the wall is taken 
into account by means of the formulation of Dorrance 
and Dore.’ Inspection of the momentum equation 
given by Eq. (4) shows that near the surface the shear 
must increase from the surface value due to the effect 
of mass injection at the surface. The increase in 
shear near the wall can be estimated by 


rT = 1,(1 + BZ) (14) 


) ro" dy + Toto (1 + B) (13) 


where Z is the velocity ratio u/u,. The shear should 
approach zero at the outer edge of the boundary layer, 
but the influence of the velocity distribution in this 
region on the skin friction is probably small. The 
success of the constant-shear assumption usually 
made for turbulent boundary layers with no blowing 
is an indication that this assumption is reasonable. 
A mixing-length theory is utilized in order to estimate 
the velocity distribution. Thus 


r = pl*|d0u/dy|(0u/dy) (15) 


Either the Prandtl assumption, / = Ay, or the von 
Karman similarity law, 1 = K (0a#/Ody)/(0°#/Oy"), can 
be used for the mixing-length variation. It is assumed 
that these formulations, originally developed for in- 
compressible flow with little x variation of fluid proper- 
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ties, apply to the present case. Both are analyzed. 
The velocity distribution is given implicitly by the 
change of variable from y to Z, at constant x, for use 


in the momentum integral; viz., 


: ol B/2 
wly” ( t °F) ia 7 
dy = (“2 Je vended | 5 : | dZ (16) 
5 Ku, pull + BZ) 


where @ is zero if the von Karman similarity law is used 
and unity if / = Ay, and 


a, = Ku, V pu Tu 
Zz 
y(Z) | [p/ pw(1 + BZ)|'/? dZ 
0 


and F and K are constants. 

The quantity a,” which appears in the above equa- 
tion is inversely proportional to the skin-friction coeffi- 
cient. At a sufficient distance from the stagnation 
point a,” should be large. Based primarily upon the 
assumption that a,” is large, the momentum integral 
with the change of variable given by Eq. (16) can be 
solved analytically by approximate methods. The 
momentum-integral equation can be put into the form 


(D ‘K*)a,?(d dx)(a,2ro"J) _ Pult-ror (1 + B) ‘Mw (17) 
where 


D =e KF 


a 
J = f fe*dy 

0 
f - [(B/ Pw)Z(1 — Z)| [B/pw(1 + BZ)\” “ia 


Successive integration of J by parts produces the 
asymptotic expansion 


dwt N n (n) c7We 

e€ (—1)" df 
J= - yy — — (18) 

Ay n=0 Ly dy 0 


For large a, with y, of order unity the dominant term 


1S 


J = (1/a,)e” (1 + B)[(pw/p.)(1 + B)]~** (19) 


In addition, if a,, is large, 


d 


(Qwe)* 7 g(x)ewe* = [(duve)2e(x)e"*] (20) 
dx 


dx 

Eqs. (19) and (20) are used to carry out integration of 
the momentum integral in the form of Eq. (17). For 
convenience, the variations of B and y, with x are 
The arbitrary constants AK and D must 
If these are assigned 


neglected 
be obtained from experiments. 
the same values as those obtained from incompressible 
experiments with no mass injection, the skin-friction 
law which results is given by 


, -= 170+ 
Cy,''? ~ 
4.15 logi [Cy Rewr(pe/pu)'~*? (1 + B)**] (21) 
where Reyr = f prettero- dX | pr 
0 
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l Dp 1/2 
= 1Z 
? J, | — a mal 


It would be fortunate if no adjustment of the constants 
were required. a” 

The effective Reynolds number Re,7 reduces to the 
ordinary Reynolds number based on wall conditions 
for a flat plate. For a conical geometry, where ry 
cx, it becomes half the ordinary Reynolds number, 
consistent with a derivation by Van Driest!* for the 
conical geometry. In all cases, including blunt bodies, 
the heat and mass fluxes, given by Eqs. (11) and (8), 
can be expected to have peak values near the sonic point 
due to the variation of the inviscid mass flux p,u,. 

The integral ¢ is discussed in the following section. 


Evaluation of the Integral ¢ 


The integrand of the integral ¢ involves the variation 
of the density with velocity ratio in the boundary 
The enthalpy variation with velocity ratio is 
Further- 


layer. 
known in terms of its boundary conditions. 
more, the pressure is constant across the boundary 
layer. For a simple gas, specification of these two 
thermodynamic properties would be sufficient to deter- 
mine the density. However, in the boundary layer 
the atomic mass fraction of each atomic species varies. 
If the number of atoms of each species could be speci- 
fied, then the molecular species that exist at a given 
enthalpy and pressure could be computed by making 
use of chemical equilibrium, and all other thermody- 
namic properties could be evaluated. 

Fortunately, with the assumption of unit Lewis and 
Prandtl numbers, the atomic mass fraction of each 
species can be determined in terms of known boundary 
conditions and a single parameter &: 


De - (Tmy = Dale + Dae (22) 
where &= (1 — Z)B/(1 + B) 


Therefore, if h, P, and & were given, the integrand of @ 
could be computed. The computation could be 
carried out numerically for each value of Z and the 
integration performed, but in view of the crude nature 
of the turbulent analysis it seems desirable to seek an 
approximate method for evaluating ¢. 

It should be recalled that the definition of enthalpy 


is given by 
1 
k=), Kf c,dT + >) KA? (23) 
i 0 i 


Suppose, for example, that the molecular species 
existing at high temperature are C, CO, O, On, N, 
and Ne, with given atomic mass fractions for O, C, 
and N. If the thermodynamic properties were com- 
puted for such a gas mixture with no condensed phase 
allowed at low temperatures, the heat of formation of 
the mixture at low temperatures would remain finite. 
In other words, oxygen and nitrogen could recombine 


but gaseous carbon monoxide would remain. A heat 
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of formation or datum for the enthalpy of the mixture 


could be defined: 


hi = >> Kph?o (24) 


where A,’ is the mass fraction of species 7 at low tem- 
perature. Hopefully, the variation of mixture en- 
thalpy with temperature and hence density will be 
similar over the high-temperature range of interest 
when the datum for each value of € is shifted to zero. 
For a combustion reaction the value of £ at which 
oxygen and fuel are in stoichiometric proportion is 


given by 
SP = B/( - B) (25) 


For graphite combustion, the mass fractions of both 
oxygen and carbon are very small at the combustion 
zone. When ~ < £*, Kc = 0. The gas on this side 
of the combustion zone, to the outer edge where é = 0, 
consists of a mixture of oxygen, nitrogen, and carbon 


monoxide. It can be shown that in this region 


al ig W = (Meo/Mc)heo"é) (26) 
When ~ > &*—i.e., 
tion zone and the surface where (Ky + Ao») & O—it 
can be shown that 
& > &* 
h® = [(1 + B*\h.” no (M co, Me) B*heo° | é = 
[B*h.° — (Mceo/Mc)B*heo| (27) 


in the region between the combus- 


The integral ¢ can be put in the form 


P\12 l D 1/2 
o> ( ) { | 5 | dZ (28) 
Pe 0 LP + BZ) 


The quantity (p/P)'/? is a function of pressure, en- 
thalpy, and & The relevant numerical data for 
graphite have been substituted for the heats of forma- 
tion and equilibrium constants. The results of calcula- 
tion for (p/P)'/? using the full set of equations are pre- 
sented in Fig. 1. The mixture heats of formation were 
subtracted out. Values of & equal to zero, which 
represents air, $* which represents the combustion zone, 
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and 2¢* which represents a fuel-rich condition were 


chosen. 
of pressure from 10? to 10~4 atmospheres. 
seen that the results fall fairly well on a universal line. 
A fit to this line was obtained by adding a correction 
term to that which results for air treated as a perfect 
gas. The expression for the fit is 


> 


te — aah)” © 0 = 
(1/H?) + 0.715 X 10-4 H¥? (29) 


where J is the mechanical equivalent of heat, y is the 
ratio of specific heats, and 


I=h-—h® (Btu/Ib) 


A slightly better fit might have been obtained with 
another form for the correction term, but the one 
chosen allows easy analytical integration of ¢. 

For graphite surface combustion—i.e., £,. < &*—and 
for other chemical systems in which a single chemical 
constituent is injected at the surface, a simple expres- 
sion for /7 in terms of the velocity ratio results: 


H = (h,, — Hg)Z + Hy. — (u?7/2)Z* (30) 


For graphite with gas-phase combustion, the enthalpy 
H can be put into a form quadratic in Z on each side 
of the combustion zone through use of Eqs. (26) and 


(27), but the slope is discontinuous at the combustion 


zone. Therefore with gas-phase combustion the inte- 
gration must be carried out separately in the two regions 
0<Z< Z* and Z* < Z < 1. Since the results in 
the case of gas-phase combustion are similar to but 
more complex than those of surface combustion, only 
the latter will be presented in detail. 

The enthalpy // of Eq. (30) can be put into the form 


A = Hyt(Z,* — Z)(Z — Z,7) (31) 

where a? = u,?/2H, (32) 
6 = (kh,/H,) + (u7/2H,) — 1 (33) 

Z, = (b/2a2)[1 + W1 + (4a2/b*)] (34) 


The superscript (+) represents the positive solution 
of Eq. (34) and the (—) refers to the negative solution. 
The integral @ can then be put into the following two 


forms depending upon the constants Z,~ and B: 
CaseI: —1/B< Z, 
Let 
sin? @ = (Z,+ — Z)/(Z,*+ — Z,7) 


k? = (Z,+ — Z,-)/|[Z,+ + (1/B)] < 1 
Then 
l ] 2 
*= UL H,aB2[(1/B) + Z,4)12 * 
- dé ‘ 
% W1 — k sin? 6 
0.715 X 10-4 X 2H,,"a(Z,+ — Z,-)? 
~ B(1/B) + ZA}? 
% sin? @ cos? 6 dé oe 
% V1 — k® sin? | - 


Computations were carried out over a range 
It can be 
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Casell: Z2,~-<—1/B 
Let 


sin? @ = (Z,+ — Z)/[Z,+ + (1/B)] 
k? = [Z,+ + (1/B)]/(Z,+ — Z,-) <1 


Then 


l 9 
re lz a BVYZ,+ — Z,)2 * 


f dé 4 
De V1 — k* sin? 6 


0.715 X 10-* X 2H,™4 : 
= ae “ [(1/B) + Z,+] X 


Ou 
(Z,* — Z,-)' f sin? @V1 — k? sin’ do (36) 
Oe 


The integrals appearing in Eqs. (35) and (36) can be 


related to standard forms of elliptic integrals as follows: 
‘ dé 

0 V1 — Rk sin? @ 

"" sin? 6 cos* 6 dé 


0 V1 — k* sin? 0 


1 2(1 — k?) 
—| — F(0, k? 
3p? | pe ( : + 


= F(6, k’) (37) 


(2 — k?) 


2 
¢ 
i) sin? 6WV1 — k*? sin? 6d6 = 
0 


1j/(1— Rr? ‘1 — 2k? 
| F(0, k?) — \ rc (0, k*?) — 


E(6, k?) — sin 6 cos 6 V1 — k? sin? i| (38) 


9 k2 


oO 


sin 0 cos 0V1 — k? sin? | (39) 


where F(6,k*) and E£(6,k*) are elliptic integrals of the 
first and second kinds, respectively, and are tabulated. 

At low enthalpies in the boundary layer the second 
terms in Eqs. (35) and (36) should be small. Then 
these expressions used in Eq. (21) would be equivalent 
to the results of Dorrance and Dore® for the skin fric- 
tion with blowing on flat plates with the mixture length 
given by | = Ay, provided the approximation of Eq. 
(19) is made in their analysis. The solution also has a 
form similar to that of Rubesin,'® who made an ap- 
proximation similar to that of Eq. (19) but took account 
of the laminar sublayer. At low enthalpies with no 
blowing the solution reduces to that of Van Driest'® 
when ¢ is inserted in Eq. (21). 


Effect of Mass Injection on Skin Friction and 
Heat Transfer 


It is often of interest to determine the influence of 
mass injection on skin friction and heat transfer as com- 
pared to that which would occur at the same free- 
stream and surface conditions without blowing. If 
the skin-friction law given by Eq. (21) evaluated with- 
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out blowing is subtracted from the same expression 
with mass injection, the result can be put in the form 


és we 
: a 1 = aln E (: + ) | (40) 
errr r 


? 


where r = Cy/Cy, B rB = m,,/(Petecl Je,/2) 
a= 1.80¢ ;,,!/? do 


The subscript 0 refers to the no-blowing case. It 
should be noted that the skin-friction coefficient for 
no-blowing depends upon which mixing-length formula- 
tion is assumed, through the term (p,/p,)'~‘”” in 
Kq. (21). 

Although large values of B can be obtained by forced 
injection, the values of 8 which result for most ablation 
materials are usually of the order of unity or less. 
Therefore, a power-series development may be more 
convenient for ablation materials than use of the solu- 
tions for @ given by Eqs. (35) and (36). The quantity 
B seems to be more popular than B for correlation of 
mass-transfer data. 

The first term on the left-hand side of Eq. (40) 
represents the ratio of two integrals. In terms of 7 and 


B they are 
l ‘ie UdZ ' 
¢ - = | = ares f UdZ (41) 
go r'/* 0 (r+ BZ)'? 0 
Differentiation of Eqs. (40) and (41) allows evaluation 
of the coefficients in a Taylor’s-series expansion for r 
in terms of B. The first few terms are 
C te (Is Ih) + ap - 
© ap 4 + See 
Cr 1 + 2a 
k '2) (Ts T;) — ag 
1 + 2a 
[(3/2) + 4a] (12o/l)? + 2a8(1 + 3a) X 
(Io/Th) + (a?B?/2)(1 + =] Xx 
(1 + 2a)3 


B?7/2+... (42) 
1 
where a -{ UZ"—'dZ 
0 


As previously noted, the integration must be carried 
out in separate regions if combustion reactions occur. 
Otherwise the integrals J, are easily evaluated; viz., 


I, = Uo = I+ 0.715 X 10-4 J; (43) 


{1 — (h,/H,,)¥? ob yh m1F 
pau i — + 0.71: 
Is (  Hy"?a? T oa f + - 
( H,'!? r b 
10 7 ~~ {1 — (&./H,)**} + = Jy aia 
oa~ adil 
a f3b — (2a* + 3b)(h,/Hy)"” 
a Hy¥!*4a4 . 
3b + 4a? 
e+ 40* 2  o715 x 
8a‘ f 
H,,'2 : 
10-4 ; 24a4 {5b — (6a? + 5b)(h,/Hy)*!?} + 
240 


5b? + 4a? 
l6a* 


ns (45) 
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where 


I l | , ( 2a* — b ) 
= : sin 
FF, “a V b? + 4a? 


b | 
sin! (46) 
Vb? + 4a’ 


| 9 / 2 
J, = —[b + (2a? — b)(h,/Hy)"?] + 
4a’ 
4, 2 5? 
c 0s i 
Sa- 


The quantities a? and b are defined in Eqs. (32) and 
(33). 

The expression given by Eq. (42) is usually adequate 
for B < 1 and therefore should be applicable to most 
ablation materials. In contrast to the laminar case, 
the relation between r and B depends somewhat upon 
both the Reynolds number and compressibility. 


Correlation With Experiments 


Much of the compressible turbulent mass-transfer 
data appear to be compared with the analysis of Rube- 
sin,!* while the present analysis is more closely related 
to that of Dorrance and Dore.’ Rubesin himself first 
compared his analytical results with those of Dorrance 
and Dore in reference 14. He found that Dorrance 
and Dore would predict a somewhat smaller value of 
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c;,/Cro for a given value of B. The difference between 
these two analyses, both of which used the Prandtl mix 
ing-length assumption, / = Ay, appears to be no greater 
than the difference caused by the choice of mixing 
length law. In other words, the analysis of Dorrance 
and Dore would tend to coincide with that of Rubesin 
if the former had been developed with the von Karman 
similarity law with / = Ay in the latter. This effect 
is due to the term (1 + B/r)*? in Eq. (40). 


Rubesin'? compared his analysis with the data of 
Mickley et al.'!° for subsonic flow and with data of 
references 16 and 17 for supersonic flow. The cor- 
relation was fairly good. Leadon and Scott'® also 
obtained good correlation with the analysis of Rubesin 
for supersonic flow. Green and Nall!’ obtained agree- 
ment with the Rubesin analysis for a transpiration 
cooled nozzle when the constant A was varied with 


mass injection in the process of reducing their data. 


For ablation studies the test facility primarily used 
is the air-stabilized arc jet. However, the power re- 
quired in such facilities to develop sufficiently large 
Reynolds numbers to obtain turbulent flow is ex- 
tremely large.”” Experiments were conducted at 
Aeronutronic’® in which graphite specimens were burned 
in a blowdown system with a controlled mixture of 
oxygen and nitrogen as the flowing gas. The Reynolds 
number per foot in the test section ranged from 1.5 X 
10‘ to 108, so that the flow was undoubtedly turbulent. 

The ablating surface was the inside of a circular cylin- 
der which formed a part of the blowdown system. 


SPEER 3474D GRAPHITE 
AT 6" FROM LIP 


THEORY «x 0.75 


——— THEORY «x 0.75 (FOR CONSTANT Cte BASED ON 
VAN DRIEST’S RESULTS FOR NO BLOWING AND 
SURFACE TEMPERATURE OF 3000°R) 


(0) EXPERIMENT P~170 psia 
EXPERIMENT P ~ 310 psia 
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CHEMICALLY ACTIVE 


The length-to-diameter ratio was sufficiently small to 
flat-plate 
method of conducting the tests 
The graphite was preheated to approximately 3,000°F 


approximate conditions. The primary 


ras briefly as follows: 


by means of an induction heater. After the heater was 
shut off the flow of an oxidizing gas over the hot 
graphite caused ignition and further heating, due to 
the heat of combustion. The equilibrium tempera- 
ture, reached after considerable time, varied from 
3,300 to 5,900°F 
parameters, pressure, and mixture ratio. 


depending upon the independent 
High pres- 
sures and high oxygen concentrations produced the 
higher temperatures. These temperatures were too low 
to cause gas-phase combustion. 

The predicted mass loss was consistently higher than 
that observed. Therefore, for the purpose of correla- 
tion the theoretical prediction was reduced by a factor 
of 0.75. The correlation with chamber pressure in the 
range of 70 to 330 psia at 100 percent oxygen is pre- 
sented in Fig. 2. Each data point represents the aver- 
age rate of surface regression for a single run at a point 
6 inches from the lip of the specimen. In this series 
of runs surface combustion occurred with B = B* = 0.75 
for 100 percent oxygen. The mass flux can be computed 
from Eq. (8). Since p,u, varies directly with the pres- 
sure, while c,, varies approximately as P~'®, the sur- 
face mass flux should depend upon a 4/5 power of the 
pressure. It can be seen from Fig. 2 that the pre- 
dicted variation with pressure is reasonably well veri- 
fied. 

Correlation with oxygen concentration over the 
range of 28 to 100 percent oxygen at pressure levels 
of approximately 170 psia to 310 psia is presented in 
Fig. 3. B* is proportional to the oxygen concentration 
I'y.. With increase in Ip, the skin-friction coefficient 
Cr is decreased due to the effects of both increased 
blowing and increased surface temperature. The 
dashed line in Fig. 3 shows the theoretical surface 
regression rate neglecting the effect of any change in 
Cy, While the solid line includes the theoretical reduc- 
tion in c,, when the von Karman similarity law is used; 
i.e., 8 = 0. Incidentally, if / = Ky (i., 8 = 1), the 
prediction would be about 15 percent higher. It ap- 
pears from Fig. 3 that the experimental data verify 
well the effect of oxygen concentration on surface 
ablation rate, including the influence of cy. 

It is not really possible to ascertain whether the factor 
0.75 is a characteristic of this particular experiment, 
or the result of uncertainty in graphite property data, 
or the result of some imperfection in the theory. How- 
ever, the results are consistent enough to warrant con- 


fidence in the turbulent analysis. 
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Thermoelastic Behavior of a Simply Supported 
Sandwich Panel Under Large T emperature 
Gradient and Edge Compression! 


CHIEH-CHIEN CHANG* 


Summary 


The treatment is concerned with the thermal elastic behavior 
of a rectangular sandwich panel which is subjected to edge com 
transverse load, and high temperature difference be 
The differential equations and boundary 


pression, 
tween the two faces. 


conditions are formulated from the principle of variational po- 


tential energy. The thermal effect on Hoff’s boundary condi 
tions is also formulated 

For the case of a simply supported panel, an analytic theory 
is presented to solve for the displacements, the deflection and the 
stresses. Some preliminary theoretical results and experimental 
verifications are shown 
an infinitely wide sandwich panel under 
cylindrical bending is also treated. For 17-7PH 
stainless steel sandwich panel with one face at 900°F and the 
other at 100°F is illustrated under both edge 


As a limiting case, 
application, 


compression and 


transverse load. The maximum stresses on both faces are given 


as a function of edge compression. 


Symbols 
A = panel area (in.*) 
ate, Birr, Cm = coefficients 
aj 
D = d/dn, operator 
(i+ f)? E’t’ ee 
Do = oes - = flexural rigidity 
1+m 1—p (Ib-in.?/in. ) 
z = Young’s modulus (lb/in.?) 
Cap = (1/2)(va.g + Us.) = strain components 
f = (t’ + t”)/2 (in.) 
G = core shear modulus (Ib/in.?) 
g = (,/G,, ratio of core shear moduli 
fl = iG, core shear stiffness (1b/in.) 
Rie = N,,,/P., edge loading parameter (positive ) 
K, = 2gR/8(1 + u) 
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(2pR/1 + w)|(1 + g) + (¢/G) 
lengths of panel in x- and_ y-directions 


(in 
al(T’ —T,) —a"(T” — T,) 
(1 + p)Do a = 
j++ TJ 
1+m ' = 
—— Do, moment of thermal differ 
‘+f 


ential (Ib-in. /in. ) 
sectional moment (lb-in. /in 
we’ 2", 


component of edge load per unit length 


face extension factor 


(1b/in. ) 
direction cosine of unit normal in a-direc- 


tion 


(i + f)? x \? t'E’ : . 
= Euler cylindrical 
1 — zp? \2b 1+ m 


) 


buckling load (lb/in. ) 
(q 2G,rRid2a ), transverse load factor 
transverse shear load in a-direction 
(1b/in. ) 
transverse load (lb/in.?) 


a(T’ — T,) —a"(T’ — T, 


(i =F 
1+ m 
P,./PVH,, core shear stiffness parameter 
(positive) 


length along panel edge (in 
thickness (in 
(i+f)/i 
reference temperature at neutral surface 
strain energy (Ib-in.) 
strain energy density function (lb/in.?) 
potential energy (lb-in 
u’/b, v’'/b, dimensionless displacements 
x- and y-components of displacement of 
outer face | (in. ) 
’/b, dimensionless deflection 
panel deflection (in.) 
deflection of the 
and ther 


w,'/b, dimensionless 
panel under uniaxial loading 
mal gradient 

a(T’ — T,) — aT" — T,) ' 

Pa (1+) = 
dimensionless maximum thermal deflec- 
tion 

52°gb/96P,, dimensionless maximum de- 
flection due to transverse loading q 

Cartesian coordinates (in. ) 

i/b, core thickness ratio or coefficient of 
thermal expansion 

b/a, panel aspect ratio 

increment of (_ ) 


determinants as specified 
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i, 7 = x/a, y/b = nondimensional coordinates 
r = sth characteristic root 
lu = Poisson’s ratio 
p = (1/2)(1+m)i 
Tap = sectional stress (Ib/in.?) 
me = shear stress of core (Ib/in.?) 
®, y, x = functions 
Superscripts 
(7) = lower or upper face 
= core 
= outer face (high temperature ) 
= inner face (low temperature) 
np = ath and pth degrees of freedom 
Subscripts 
a or B = free indices each taking x or y indicating 


direction 
,a or ,B = free indices taking x or y indicating deriva 
tive with respect to x or y, respectively 
f = face 
= ith direction 
derivative with respect to x or y 


ll 


x or ,y 
L = loading 
= temperature 


Signs of Quantities 
Mas, Nas, Qa, 7 are positive as shown in Fig. A-1b 
U is positive when work is done by the panel on the boundary 


restraint. 
lV’, the potential, is negative when work is done by the positive 


external loadings on the panel. 


(I) Introduction and Statement of the Problem 


(1.1) General Sketch 

INCE THE PIONEER WORK on honeycomb sandwich 
S panels developed by the Glenn L. Martin Com- 
pany during the World War II, sandwich structure 
has been widely applied to high-speed aircraft and 
missiles, and will be increasingly applied to the struc- 
tures of future re-entry and space vehicles. Its high 
strength-weight ratio, its smooth and firm external 
surface, and its good heat-resisting and insulating 
properties are some of the promising features which 
will be desirable in future applications. 

One of the important types of sandwich structure is 
the rectangular flat sandwich panel under edge com- 
pression and exposed externally to very high tempera- 
tures because of aerodynamic heating at hypersonic 
speeds. Under such high temperature differences be- 
tween the inner and outer faces of the panel, thermal 
stress and warpage become challenging problems of the 
space age. The present paper will bring some analyti- 
cal answers to these problems. 

In references 1 and 2 have been presented treatments 
of the simply-supported rectangular panel under edge 
compression, lateral loading, and initial curvature, with 
the core orthotropic and the faces of different thickness, 
Indirectly, the earlier 


material, and temperature. 
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theory can give approximate solutions for the sandwich 
panel with a high temperature gradient across the 
thickness provided that the temperature in each face is 
uniform. The present theory takes the thermoelastic 
effects directly into both the differential equations 
and the boundary conditions, and may be considered 
as an improvement of the early works. For the sake 
of brevity in the present paper, references 1 and 2 are 
quoted frequently. 

The method of variational 
treated in the literature, as shown in the references. 
Appendix A gives a short resume of the variational 
principle as it is applied to the thermal stress problem 
of sandwich panels, including all possible types of 


energy has been well 


boundary conditions. 

The present treatment is mainly for simply supported 
sandwich panel with the core orthotropic and the two 
faces different in thickness, material and temperature. 
It is found that one nonhomogeneous term occurs in 
the differential equation in the direction of deflection. 
There are four nonhomogeneous boundary conditions. 
The present approach is the first to show that the solu- 
tion for displacements is separable and that there are 
an infinite number of eigenvalues or degrees of freedom. 
Thus, the general solutions can be expressed as infinite 
series. Trigonometric functions are chosen, in general, 
to satisfy the twelve boundary conditions for each 
The general method and solutions 
However, 


degree of freedom. 
of the displacements and stresses are given. 
the theory is much more involved than the previous 
treatment of references 1 and 2. Some preliminary 
theoretical results and experimental verifications are 
also given. The numerical computation is now in 
process. The complete results will be reported in later 
papers. 

The cylindrical bending of a sandwich panel with 
thermoelastic effects is treated as a special case. The 
results are made equivalent to the general solution by 
letting the unloaded edges be infinitely long. An ex- 
ample is given for an infinitely wide stainless steel 
sandwich panel with the outer hot face at a tempera- 
ture of 900°F, the inner cold face at 100°F, and a 


compression load along two opposite edges. 


(1.2) Assumptions 

The basic assumptions are: 

(a) Two faces are isotropic but may be different in 
material, thickness, and temperature. The Poisson's 
ratios of the faces are the same. Bending rigidities 
of the individual faces may be neglected. 

(b) The core, being weak, can be orthotropic in trans- 
verse shear rigidity, but approximately homogeneous. 
The average temperature effect on the core shear 
modulus can be considered. 

(c) The bonds between the faces and core are strong 
and no bond failure is considered. 

(d) Uniform edge compression is applied on opposite 
sides of the panel at the neutral plane. Transverse 


loading is light. 
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(e) Temperatures at the two faces are different but 
constant. 
(f) Small deflection elastic theory is considered. 
Justification of the above assumptions is given in 


reference 1. 


(II) Analysis 
For convenience, the analysis on thermoelastic be- 
havior of a simply supported sandwich panel under 
loadings is described in the following steps. 


(2.1) Differential Equations of the Problem 
The general derivations of the differential equations 
and boundary conditions are given in Appendix A. Fig. 
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A-la shows the simply supported sandwich panel under 
transverse loading, edge compressions and shear. Fig 
A-1lb shows the loadings on an element of the panel 
The general differential equations and boundary condi 
tions for a loaded sandwich panel with thermoelastic 
energy 


effects are obtained from the variational 


method. 
The following assumptions are imposed : 


(i) The panel is simply supported. 

(11) Edge loadings N,, = N,, = 0. Now, N,, is 
in compression and so is negative. Introduce 
Nw = —N,,. Then N,, > O when N,, is in 
compression. 

(iii) The outer face temperature 7” 
than the inner face temperature 7”. 

(iv) The transverse loading g is constant. 


is much larger 


Corresponding to the above assumptions, the system of three differential equations of (A-11) can be reduced to 


L (1 + m)t'E’ 


12" 2s + (1 = pat” wy ote (1 + w)v" ry] 


l = l : 
+ = A. ( a so u’ + fw,’ = () 


(a) z 
2 (1 — pg?) t 
_Llo+m’Eé , : ' l+m,/fl+-m,.,, 
(b) (20 yy t+ — we er + +4)’ ,y] — ; f. i v’ + tw, =O (1) 


2 (1 — p?) 


Pm a a, — Pee eee a ae ai; aa 
(¢) -ift, ( ; uw’, + tw m — tH, ( 7 Pw + tw » + Ny’ wy = 4 





where 
(a) (= (4+ /S)/,, [f = + t")/2] (b) m=tE'/t"E"| (2) 
(c) A, = iG,. H, = iG, f ‘ 
The equations in (1) can be further reduced to nondimensional forms. Define 
(a) £=-</a, 7 = y/b (3 
, , / 0) 
(b) u’/b = 4, v'/b = 2, w’/b = wif 
With the introduction of Eq. (3) into (1), the corresponding system of nondimensional equations is: 
(a)  (4-e/m*) [287 cg + (1 — ptm + BCL + #)2%9] — 2gu — (BgPa/p)w.. = 0) 
(b) (4r-/m?) (20,4, + B21 — uw)? + BU + w)u,,] — 20 — (Pa/p)w,, = 0 (4) 
(c) —2pgBu, — 2pv,, — gP?B2aws, + (rk — l)al*w,,, = g/Gy \ 
where 
(a) a =i/b fe) P= G@+/f)? t’E’ (=) 
(b) B = b/a 1— wp? 1+ m \2b (5) 
(ec) g =G,/G, (f) Rie = Ny/P. | 
(d) = (1/2)(1 + m)i (g) +. = P./@H, 


The first four nondimensional parameters are understandable as defined in the nomenclature. 
As discussed in reference 1, it can be shown without much difficulty that ?, 


of P, in (5e) needs some explanation. 


is the Euler load per unit width for cylindrical bending of the panel as the width goes to infinity. 
measure for the applied edge compression N,, and for core shear stiffness /7, = 1G). 


The introduction 


It is a meaningful 
Thus, the nondimensional 


parameters k;, and r, are introduced respectively as shown in (5f) and (5g). 


(2.2) Boundary Conditions and Hoff’s Modification 


For the case of simply supported panel, the boundary conditions in (A-12) can be rewritten in nondimensional 


forms as follows: 
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(a) Bu, + wv, = R é= + ! 
(b) v,, + Buu, = R "= *1 
(c) uw, + bv, = 0 &é&= += 1 
(d) uw, + Bo; = 0 7»= *1 (6) 
(e) w= 0 gé= += ] 
(f) w=0 7= +1 

a’(T’ — T,) - a’(T” — T,) : : ; 

(g) R= (1 + wp) is a nondimensional constant 


1+m 


It is difficult to satisfy exactly all the boundary conditions. In reference 6, where the sandwich panel was loaded 
at uniform room temperature, Hoff showed that (6c) and (6d) can be relaxed to 


v=0,& = Fl; w 0, 7 + | 


However, for thermal stresses Hoff’s boundary conditions should be replaced by 
R 7 
v= n (é # 1) (6°) 
Us 
e 3 l (6d’ 
u= E (n > 4 ) x ) 
B(1 + pw) 


as Shown in Appendix A. 


(2.3) Choice of General Forms of Solutions 


Initial choice of the general form of solution to satisfy all the boundary conditions can save a great deal of analysis. 
3y inspection, the reasonable choice of solutions is of the form: 


R . eye, pr 
(a) u(&, n) &é+ > > A, sin — cos n 
B(1l + pw) n=1 p=1 2 2 
R ct nx. |. pr 
(b) v(& n) = n + >> > A»,"” cos — sin n (7) 
] + M n 1 p=1 2 2 


- nT pr 
(c) w(t, 7) = > yy A;”” cos r* £ cos i 
n p 


where » and p are the eigenvalues or degrees of freedom. All the boundary conditions are fully satisfied by the above 
equations provided only that the odd values of 7 and # are taken. 

Some insight into the behavior of the dimensionless displacements wu, v and the dimensionless deflection w (here- 
after simply referred to as the displacements u, v and the deflection w) may be obtained from the nature of the load- 
ing. It can be seen that u(é,7) must be an even function of 7 and an odd function of —. Similarly, the displacement 
v(£,9) must be an even function of € and an odd function of yn. Finally, w(,y) must be even with respect to both & 
and y. Certainly, the choice of (7) satisfies these requirements. 

Substitution of (7) into (4) yields the following three expressions: 


ip ] np . n y * 
DD (bli? + xd” + ids”) sin — Ecos 9 = Kik 
n=1 p=1 y : 


am a os nt _ pr , 
> 2 > (Po \" + xK2A 2 + PoAs *) COS > é sin : 7 = Kon (S) 
n=1 p=1 Z ys 
nT pr 


> ie (Pz. | ‘a + X3. 1,"? + 3. i" )cos : © cos ‘ n K; 
n=1 p=1 7 7 


where 
(a) @ = —7,(1 — p)p? — 2677.n* — 2p (g) bd; = —Bpgzrn 
(b) x1 = — BU + 4)npr, (h) x; = —prp 
(c) V1 = (Brgl?an) 2p (1) v3 — (al2r? 4) [B2gn? _ (rR, = 1)p?] 
(d) & = —6r.(1 + w)np G) Ay = 2¢R B(1l + p) (9) 
(e) x2 = —2r.p? — B'r.(1 — un’? — 2 (k) Ke = 2R/(1 + uw) 
(f) Yo = (Pamrp)/2 - 2oR ( 
vs p)/<0 (1) As = - at+eat+s 


l +e Gy, 





484 JOURNAL OF THE AEROSPACE SCIENCES—JUNE 1961 


(2.4) Expansion of Nonhomogeneous Terms in Eqs. (8) Into Desirable Double Fourier Series 
The nonhomogeneous terms in Eqs. (8) should be expressed in terms of the corresponding double Fourier com 
ponents which appear on the left-hand side of each equation. This can be done by using the classical method of 


Fourier analysis. In short, 


(a) Kyé 


re pr 
p a Be By” sin — cos n 
a ; ; 


2 2 
2 , nw . oe 
(b) Kon uu > B,"” cos - & sin a (10) 
i p ~_ _ 
. in b un 2 1 
(c) Ks } a > B;”” cos 5 § cos 2 n 
n=1 p=1 zZ 2 
where the Fourier coefficients are 
(a) By" Zp Ber - j 0) (n, p even) 
3n ((32K,/2r'n2p) sin (nw/2) sin (pr/2) (n, p odd) 
t ] 
(b) Br = pe =) OO . iden (11) 
; ((32K2/r*np”) sin (nw/2) sin (pr/2) (n, p odd) 
I 
j () (1, p even) 


cr + PPP 4 ype 


(c) B;"” 


((16AK3/m°np) sin (nw/2) sin (pr/2)  (n, p odd) 

provided gis a constant. Here 4”” is introduced for B,”” in each of the three coefficients so that the notation can be 
simplified in later treatment. For the case in which g is a function symmetric with respect to — and n, 5;"” can also 
be obtained without much difficulty. 

(2.5) Evaluation of Coefficients 


If Eqs. (10) are substituted into (8), we have 


. UT Da ) 
(a) } 8  # (BA 1"? + x1A2"” + piA;"? — B,"") sin > § cos Pe oy 
n=1 p=l1 2 é | 
: a tls nx. . pr 

(b) do do (42.Ay"? + x2A0"” + WoAs”? — Bo”) cos > § sin e y= 0 > (12) 
n=1 p=1 = “ | 
nT. b ar | 
(c) >. 23 (dA 1"? + xX3A o"? + W3A 3"? neal B3"”) cos > § COs n= OV 
n=1 p=1 “ i ) 


which are true for all values of £ and y from (—1) to (+1). This means that for each degree of freedom, p and n, the 
terms in each square bracket must be identically equal to zero. In other words 


BA yj"? + x1A2"? + YyA3"? = By" | 
B.A ,"” + xX2A2"? + PAs"? = B,"” (13) 
3A "” + x3A2"? + YxA3"? = Bs" 


which is a system of three simultaneous, linear, nonhomogeneous algebraic equations. Rewriting (13) in matrix form 





gives 
; X1 Yi A," By"? 
P, X2 Yo A 2"? | = B,"” ( 13’) 
®: X3 V3 A 3"? B3”” 
By means of conventional techniques of matrix algebra such as are given in reference 12 it can be shown that 
A"?A"? = A” ~~ (i. = 1, 2, 8) (14) 


where A”’ and A,"” represent determinants such as 
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P; x1 V1 alta? E n282 n2B2\2 
(a) A”? = |, X2 ¥2| = pr. Se fee (eed 2 oe r.p* + 
. vs . al pe 

2 n*B* i aby ,. ‘ j ts 

2p? | g + - — 2p? (¢-—1)? —k& 4) 1 —w) {1 t+ r.p* + 
p- p ( p- 

f., ns? n?3? . } 

2p* 1 2 + = ret] (1 — pw) (1 + 2g r.p* + 27 
p p / 


; np s 
x v1 iy al*x® p P ne? . 
(b) A,” = | x2 » Br) = B» -| 2p4r.2ky. | z — — 2p4gr, (1 — ra 
1 Xx t g 9 f 5 

we y: BP +t np Pp’ 


=e 


, Pang 
B°n*p*r.*kie(l — w)(1 + 2) + 2ptr bs + c= xX 


2p 
E B?n? 
Ea = wy (: ae °) — etl — £) 4 28 
p* 
(15) 
eb) "Dp * 
P, v1 By al? ae 2 B?n? F Bin4 
(Cc) A"? = Py Yo Ba”? = a >= 2p'r.*ki, _ “ T 2p*zr, l1— : + 
®, Ws B,"” f p- p 
Pap 
pr *ky (1 i u)(1 + 2) + 2p | wm cn a a 
zp 
- B°n? a 
Ea — p) (: + ¢g re ) + 28*n’r.(1 — g) + 28 
2 
c ‘ > np 
BD, X1 ~ ies : : B°n? 2 i: Bn? 
(d) A;”"? = | @ X2 3a"? | = B"?prp | 7,2p*(1 + gZ) (1 — w) (1 + = + 4gp*r,{ 1+ 
Pp, x3 B,” p p- p- 
: z pen sn" \* 7 _ B’n* B7n? 
Ont orl = 99" 1 1 > + 2r.p*(l — w) {1 + €—, ] + 47. (2 + 2 ae 
| p- p- 
g 1% . av .. te 
(e) C” = 2 ———~ (oie ia’ 
G, «np 2 2 
From matrix algebra 
A”? = A,*?/A” (¢ = 1, 2, 3) (16) 


Eqs. (7) with the coefficients (16) are the s« lutions of the problem. 

Full development of the theory will be given in a later paper. Some simple examples may be of help to the reader 
in estimating the maximum deflection of the panel as a function of loading. Let it be assumed that the transverse 
loading g = 0. The maximum deflection can be evaluated from (7c) if A;”” is calculated from (16)—.e., 


A3"? = A3"?/A" (17a) 


Now from (15a) and (15d) 


is afte? - (: na n?3? (1-2) (: * n2B2\2 rp) + 2p? (< 4 a) — ( 2 - 
= , pr. g p? hb) p? ef 2p g p? 2p > (g — ) — 
n*3°\? n*B? 232 
hi. {| — n) (1 + =) rep! + 2p? (2 + re )| re + t — x) (: +3 a ) rp? + 22 | (17b) 


64Rp p nx . pr | ( B?n? ( B?n? 
‘ A;*? = r, sin sin 46+yr, p214+ 2)1 -— 1+ Pd 7 
ind “ll +a) e S 9 > +r p g M) p? a" ) (17c) 


To determine the critical value of k;,, A”” in (17b) must be set equal to zero, which yields 


(: ea at (; 4 or) es he ai (« 4 | ‘cos (“) — 
ag _ ) Dp* 7 — 2p? =» 2) 
ave Ku p cb b : p? p p ( g) 


ky, = ( 17d) 


n*B?\? St. . oF _ np? ; 
t w(t + *) rp* + 2p (z+ )| V. +{a -w (1 + @ r) rat + 28 | 
p* p* Pp 
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Fic. 2. Comparison of theoretical and experimental center 


deflection distributions w,/i#, along the line & = O at different 
values of loading parameter kj.. 


Therefore for each particular panel, A;”” can be calcu- 
lated if the detailed specifications are given. For in- 
stance, Table 1 gives the specifications of the chemically 
bonded panel which was donated by Boeing Aircraft 
Company. With this detailed information the values 
of R, %, A;'', As)8, A;*!, A;*%, ete., were calculated and 
are listed in Table 1. Just for demonstration, only 
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approximate calculations were made. In Fig. 1, the 
dotted curve is from experimental data and is plotted 
nondimensionally with k;, = N,,/P, as loading param- 
eter and w,/w%, as deflection parameter. Here w, is 
the maximum deflection of the panel and @, is the maxi- 
mum thermal deflection for cylindrical bending. In 
fact, 
wm, = (1+ m)R/2ai 


can be calculated and is shown in Table 1. The solid 
curve was calculated from (7c) up tom = 3 and p = 3. 
The comparison is satisfactory for kh, < 1.3. For 
higher k;, the panel deflection is too large for com- 
parison with the small deflection theory in this paper. 
Besides, plastic effects may be so predominant that 
the present elastic theory becomes questionable. 
From the test data of reference 16, the deflection dis- 
tribution (at & = 0) along yn at different loadings is 
shown in Fig. 2. The deflection distributions at cor- 
responding k,, were also calculated with the present 
theory. The comparison of the experimental and 
theoretical values of local deflections is satisfactory 
as long as ky < 1.077 (corresponding to N,, <500 
Ib/in.). For higher values of k,,, the comparison is poor. 


(2.6) Evaluation of the Local Bending Moments, 
Transverse Shear, and Stresses 


Knowing the solutions for u, v, and w in (7), the local 
bending moments M/,,, \/,,, M/,,, and transverse shears 
Q, and Q, are given by 


(+ fyE't’ 


M,. = —— [(Bu, + wv,) — RI 
| — S ’ 
(+ fyk't’ 
M,, —= ; es [(uBu ¢ + v ao —_ R| 
(+ f)k't’ 
M,, = : (4, + Bu ¢) (1S) 
2(1 + p) 
=. ; 1+m 
QO, = Gi (sie + u) 
Qa 
a 1+m 
QO, =Gl (1, + :) 
Qa 


The local stresses are obtained from the equations 


Ter = M,,/(i + f)t' 

Ten [M,,/(@ + f)t’| + [N,,/t'A + m)| 

tr, = M,,/(i+ ft’ (19) 
?., = O,/t 

7, = Q,/t 


Thus, the method of analysis is completed. 


(III) Cylindrical Bending Case and an Example 


The complexity of the above analysis justifies ex 
amining the simplified case of cylindrical bending 
i.e., the loaded sandwich panel becomes infinitely wide, 
ora — o. In this case, « = 0 is assumed. This 
means that the supports of the panel on the loaded 
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Fic. 3. Ratio of maximum deflection to thermal deflection (w,/t,)n-» of cylindrical bending vs. k, at different values of r- 
edges are so rigid along the loaded edges that no lateral at 7 = + 1, and (4’) is reduced to 
displacement u is allowed. Thus, v and w are inde- A(D)v = 0 A(D)w ~9/6 (20) 
pendent of x. Also,8@—~ 0asa— o. The equations 
in (4) degenerate to d 
: where D = and 
P 9 dn 
(a) (47,/m?)v,., — v — (f/2p)aw, = 0 (4’) 
= j2 es ho 7 ( : “or 
(b) 2pv,, + Park, 1)Win = 9 G,/ A(D) = —(4r,f2/x?)aD?((1 — r,k1,)D? + 
where = 1 is chosen and g is either constant or a (m?/4)Rie] (21) 
function of y. The boundary conditions in (6) are ie nifittoen of cements tims enbe 
reduced to 
. > 7, R+ Pal(Q/p) . PaQ 
(a) v, = RK (bt) w= 0 (6°) (a) v= ¥/P! sinh Ain — n 
Ai cosh Ay p 
i , R + Pa(Q/p) » 
TABLE | (b) w= vp cosh Ayn — (22) 
Specifications and Calculations for Boeing Sandwich Panel aM cosh dy 
; t?a 
Material Aluminum alloy o( + ') — (R/a,r;) + Or? 
2024-T3 apa; 
fa f 0.008 in. 
i 0.375 in. where 
2b = 2a (square) 12.0 in. 
E’ = FE” 10.5 & 108 psi . —— 9 roma? 2 ) 
T’ (average 245°F (c) aq = (ad: /2p)/[1 (4r,A,°/ 4 | 
T” (average) 80°F . ° 2 . 
T,, (average) 162°F , (d) A? = —W*k,./4C1 — rRie) or ; (233) 
G, 54,000 psi Ai = ir,’ = (im 2) (Ri. (l — rRie) |' : 
Gy 32,400 psi 
gq 0.6 os 5) yy 72 
a’ = a” 13.5 X 1076 in./in.°F (ce) Q = 9/(2G,rhid*a) 
(Eq. 5g) 0.022 . , : : 
P, (Eq. 5e) 464 Ib./in. Besides the edge compression N,, acting on the 


R (Eq. 6g) 

w, (Eq. 24d) 

A;!! (Eqs. 7, 17a) 
A;!3 (Eqs. 7, 17a) 
A;°! (Eqs. 7, 17a) 
A;33 (Eqs. 7, 17a) 
(Rie)er (Eq. 17d) 


1.45 & 1C€ 

0.0227 

(0.988 — 0.262k;.)~ 
—(14.4 — 0.296k:.)~ 
—(15.3 — 1.67 ki.) 

(39.5 — 3.40 ki.) 

3.78 





panel, there are two more types of loading which con- 
tribute to the displacements v and w—i.e., KX and g, 
the thermal and transverse loads. Sometimes, it is 
convenient to decompose w into two parts, 
w, due to R and w». due to g—1.e., 


namely, 
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(a) w= wy + w ) edge compression, w;,/%, in (24b) is plotted for the 
middle of the span (y = () in Fig. 3. 
where 


. | 8 1 feos ry’n (c) 
(b) w, = a| =< ( ~ - 1) | . 
wr” ky. \ cos \y 200 1 
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f (24) | ao ss | 
1 (= di’n | 160/—_ 3 effects Nyy. 9274 «T= «'T>0 
-1)-G-=n|, a 
ki. \ cos Ay’ | a 
} 140e— ---- 2 effects Nyy ond K'T~&"T'>0 
. 
(d) my = 
* 
a’(T’ — T,) — a" (T” — T,) , S120 ++} — 1 ‘ 
one (1 + w) | . @? 
2ta ° ] | | | | 
° Ler 
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(e) tw. = 5mr7qb/96P, ré |f_| | ° 
: ~ ar | 
Eq. (24b) shows that the deflection w, at any fixed 80 | Z | 4 
position along the span is always proportional to w%, t t | | ett Tis 
° ‘ = ~ r . “ Y | | 
which is given in (24d), and is the maximum deflection 60 | x4 +--+ a ‘ 
. ® 4 | vy - 
per unit span due to the temperature difference of the ee 1} — +1 Le4 
| - | 
faces. Inside the bracket, k;,, which is defined in (5f), 40 |“ iss 4 “al 
- : : : ~; af 2 
is the ratio of the applied edge compresssion JV, to the A Lae" 
“ 
Euler buckling load P, as given in (5e). The char- 20 v7 ee jan 
acteristic root Z ~=" | 
a + } ee + 
= /9 / = 1/2 ° a aes I : 
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is a function of the core shear stiffness parameter, 7,, — : 
: é Fic. 5. Maximum stresses of stainless steel sandwich panel under 
and k;,. For a given face temperature difference and loadings, N,,, q and temperature difference. 
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For a perfect rigid core (r, = 0), w;/%, approaches 1 -Nyy Nxy 
when k;, — 0-—-i.e., no edge compression. As hk, in- == ee bE ; + | 
creases from zero, w,/%, increases from unity and is + “Nyy —H A 
shown as the first curve on the right in Fig. 3. As “Hoe ~x 2b “4 ag 
soon as k;, approaches unity, w;/% approaches infinity. nas ' Nyy |} 
That is, critical buckling failure occurs if the panel is _ is Ey ‘xx LL Ae 
still elastic. This curve applies not only to cylindrical » ee gg 4 : 
bending of a sandwich panel of perfect core shear i bodbduddddd dd) 4 
stiffness but also to cylindrical bending of a plate. lt siaasiaaiaiandietataecmena Xx 
This curve checks with the result in Eq. (27) of Chapter t' fr 
] in reference 14, although the physical interpretations MY 
are different in the two cases. 

For core of imperfect shear stiffness— i.e., 7, > 0O—a 
family of curves are similarly plotted as shown in Fig. 3. 
The critical buckling load (N,,)-, or (Rie)cr decreases ——— Che ae Te 


as r, increases. 

As far as the contribution of transverse loading qg to 
the deflection is concerned, (24c) shows the analytic 
expression of w. which is plotted in Fig. 4. It should 
be noted that # is the maximum bending deflection 
due to a uniformly distributed load g. For the case 
r, = O, W. is given in Eq. (15) of Chapter 1 of reference 
14. 

The most useful expression is the bending moment 
M which is related to the derivative of v as follows: 


(a) My = [E’’/(1 — w*)|@+f/)@, — R) 


4gb? (‘= Ai’n ) 
Mo = 1 
( * | cos ),’ | 


[E't’'/(1 — w2)]@ +f), — R) 


4 qb? (‘ cos \y’n ) ‘ 
= M, ” — ] f 2: 
( + =) cos J,’ oan 
(c) My = (1 +4) X 
a'(T’ — T) a a"(r" =a T,) a 
t+ f 
(d) Do = [(@4+f/)?/0 + m))[E'’/(1 — w?)] 
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Fic. 6. Shear stress 7y: of the sandwich core vs. Rass 
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Fic. 7a. Yop: Simply supported sandwich panel under edge 
compressions (N;z, Ny,), shear (N-z,), and transverse loading (q) 

Fic. 7b. Bottom: Resultant forces (Nzz, Nyy, Ney, Qs, Qy) 
and moments (M,,, M,,, M,z,) acting on a sandwich element of 
unit area. 


The transverse shear can be shown to be 


O 


eV 


= —(2/m)*(77,/a) [(1 + m)r,R + 
(q/G,A)(1/Rie) |Ax'(sin Ay'n/cos Ay’) (26) 


In (25c) M, is the moment of thermal differential 
and may be interpreted as the moment required at the 
panel edge to keep the panel flat when subjected to the 
face temperature difference, provided both V,, and gq 
are zero. If g remains zero, the coupling between Mo 
and NV,, is given by the terms in the right-hand paren- 
theses of (25a). The remaining part in the equation 
shows the coupling between finite values of g and N,,. 

The important stresses can be calculated by substi- 
tuting (25) and (26) in the following equations: 


@) 2) = Mo = 
a) te  104+f)' #0+m)” 

M, m " (27) 
b) ty’ = — ——"— + — K 
(b) tw t”(? + f) + u“(itm) ™ 


(c) tz = Q,/t 


Illustrative Example 

As an application to the technical problem, the fol- 
lowing example is shown. Consider a sandwich pane 
made of 17-7PH stainless steel face sheets and square- 
celled core with the following structural details: 
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t' = t” = 1/80 in. T’ = 900°F 

i = 0.500 in. T” = 100°F 

2b = 12 in. Tl, = 0 F 

E’ = E” = 30 X 108 psi q = 100 lb/ft.’ 
G = 50 X 108 psi P, = 3,409 Ib/in. 
a’ = a” = 6.5 X 10-* (°F)— t, = 0.44122 
a,’ = 100 X 10% psi (yield stress of hot face) 

o,"= —180 X 10% psi (yield stress of cold face) 


The stresses 7,,’ and 7,," vs. edge compression N,, are 
calculated from (27) and (25). It is found that for the 
present structure, the cold inner face is more critical 
than the hot outer face, as shown in Fig. 5. The ulti- 
mate load should be slightly above 1,500 ib/in., at 
which point the yielding stress of the cold face is 
reached. The stress r,,’ of the hot outer face is still 
below the proportional limit at N,, = 1,500 Ib/in. 

The deflection w corresponding to the ultimate load 
can be calculated from (24), but the accuracy may not 
be satisfactory because the panel for N,, = 1,500 Ib/in. 
is far beyond the elastic limit. Further improvement 
on the estimation of maximum panel deflection will be 
discussed in a later report. The stress 7,, vs. edge 
compression N,, is calculated at 7 = 1 and plotted in 
Fig. 6. 


Appendix A 





The Duhamel-Neumann thermal stress-strain relations 


for plane stress are 


l+u m 
Cap _ E@ ig Tag” — 1 a Ty _* hat + 
¢ m 


a? (TO — T,)be, (A-1) 


It is shown in the small displacement theory of elas- 
ticity (reference 13, Article 62) that if a strain-energy 
density function LU) exists, then 


dUo = Tap deag 
and, consequently, 

Tap = OU/0egg (A-2) 
Inverting Eq. (A-1) we find 


EO m 
op” ee nies a (i) +a 2a a D5, = 
7” a (eu "t=." :) 
(TO — T; EE 
ee a 
~~, 
From Eqs. (A-2) and (A-3) 


y 7 (i) 
OUo _ hi — fe) + a. - ¢. (05 os 
Cas? l+u ap i-, ” ap 

ar 





a(t — THE” 


etal , , , 1— 7 
Derivation of General Differential Equations “ 
and Boundary Conditions of Sandwich Panels; Integration of these differential equations yields 
Principle of Minimum Potential Energy 
7 (*) m ; 
Energy Integrals UO, = — (cas) + —— eyes) Cag’? — 
myles ; 2(1 + un) is 
Figs. 7 show the loaded sandwich panel and an ele- . 
ment of sandwich panel under a temperature gradient. a ( pec — T,)E® e, (A-4) 
_ ee i bu TT 
The strain energy of the faces is 
; E® im a%(T® — T,)E . 
VU, = (© ff UjdA = 0 ff} ——~ | ¢ag + —— e,, 856) €ag? — ————_———— (or dA 
J A Pye 4. u) B 1 - YY 8 3 1 — yp YY 
(A-5a) 
where Cap? = (1/2)(uag + Ugo”) (A-5b) 
Here conventional notations of Cartesian tensors are adopted, with the understanding that u, 4” = (0/08)u,“” is 


an example. 


The strain energy ‘n the core can be written as 


Each symbol is defined at the beginning of this paper and no further explanations w-Il be given. The 


contribution of bending energy in both faces is neglected. 


. i . f, 1+m : 
U=- f ¥ 4G (1. —— ue! dA (A-6) 
— A a 
where m= t'E'/t"E” = —u"/u' = —v"/o' (A-7) 


is introduced under the assumption of equal Poisson’s ratios for the two faces. 
Finally, the potential due to the edge load NV, and the transverse load g can be shown to be 


V=- ff (¥., a + ru) dA (A-8) 
A 2 


According to the principle of virtual displacement, a necessary and sufficient condition of equilibrium of any elastic 
body is that the work done by the external loads (— V), less the changes in strain energy (U), must vanish for every 


Variational Principle 
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virtual displacement—i.e., 6(U + V) = 0. In other words, summation of all the increments of the above energy 


integrals (A-5), (A-6), and (A-8) due to any virtual displacements must be identically zero—i.e., 


ov F ve ( 
0 = (U + V) = 8U/ +8U,’ + 80 + 6V -ff f +s (capeas +—= euie) _ 
l+u Sed | 

AP ~T)< wt? ~ Td. , 
Pet Be | dA +1 ffx ye G. (ta, eee ue! 5 (' CF ee ue! i 

=. A t . t 

{J (N as q/dw’ + giw’) dA (A-9) 
A 


where the first integral on the right-hand side is obtained by means of (A-7). Using the divergence theorem, it can 


be further shown that (A-9) becomes 


(1 t'E’ rz’ 
0 =U + V) =f{ fi- + m) (c ase +3 Hes -) +; + ~ [a'(I" — Ty) — a"(T" - T,)|qbéu¢'dA + 
\ 1+u ~~ - f 
(1 Hf, {\ / ’ = l . 
{Jz us ” “( s = u', + tw’ \e A ff i, > A. ( s = teat iw’ aa — 


(1 + m)t'E’ ; 
q + Nast aot b6w'dA + f J sms + - ex) ne | — 
= 


s ( lth 
a'(T’ — T,) — a"(T” — T,) ,,, - |l+m Fr —— . 
- t E'n.$ bu’, ds — # IS tl. = ua + tw’ a | Ma t+ Nagw ink 6bw’'ds (A-10) 
l — w f S { a t : : f 
where s is the perimeter of the area A, 7, = G.i, 1, = G,i, and n,, ng are components of the unit normal to the 


boundary of the area A. 
(A-10) is true for any arbitrary changes 6u,,’ and 6w’. 
Thus, three differential equations are obtained respectively in the direction of x, y, and s: 


This means that the integrand in each big bracket must be 


identically zero. 


—  —(1 + m)t'E’ Z l+m,/fl+m , _ 
(a) : [Qt er’ + (1 — wy’ + CL + w)o,2y/] + i i. ( ; u’ + iw,’ = | 


2(1 — mu’) 
Juz = aad = 
— [a’(T" — T,) — a (T”" — T,)}.2 
i—<~¢ 
—(1 + m)t'E’ ; Less li+m,/f/l+m,., 
(b) = [Qe yy) + (1 — weer’ + (1 + w)t zy | 4 - i, ( - ve + iy! (A-11) 
2(1 — pu?) : t l 
‘rt 


_ [a’(T’ — Tn) — a’(T” — Tn)\, 
L=- 


-{f{/l+m : se lm j = , os . 
(c) —tH, ( — a fw." = i, ( j v mY +f W yy ‘) + N,2W, rz yo ZN x80 cy + NA yy yy q | 


The above are identical to Eqs. (4) 


where the subscripts ,x and ,y refer to differentiation with respect to x and y 
The complication of the thermoelastic 


of reference 1 except for the nonhomogeneous terms on the right-hand side. 

sandwich panel is due to these nonhomogeneous terms as seen in the treatment. 
The line integrals in (A-10) give the boundary condi- 

tions and cover all the possibilities for a loaded panel. 


Boundary conditions for a simply supported panel are: These conditions are sufficient to solve the problem 
(a) te ’ if g, T’, and 7” are given functions of x and y. 
a uy Vv, = , p % ; ‘ 
gs - a eee enna It is convenient to write the displacements w’, 
(T’ — T,) — a®(T” — T,) - 
a (1 + yp) v as 
m 
tor x= £¢ u = up) (x) + uz” (x, y)t (A-13) 
*) (2) Gi A-15 
(b) vy’ + uu, v = o7‘”(y) vz,‘ (x, y)f 
'’(T’ — T,) — a’ (T” — T,) ‘ wae hea ~ 
n SF i+ | (A-12) where the subscripts 7’ refer to the thermal origin and 
i = sic subscripts Z refer to the mechanical origin. Since 
y = + b Tr 7 
' up (x) = aM (T® — T, |x 
(c) uw,’ +2,’ = 0 x= +a . | 
(d) uy’ tv,’ = 0 y= +b) therefore 
(fe) w’ = 0 x=+ta ur'(x) — up”(x) = (1 + m)ur'(x) X 
(f) w’ = 0 y= +b fa’(T’ — T,) — a"(T” — T,){ x 
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where Eq. (A-7) is used. A similar expression can be 
written for v7. Rewriting Eq. (A-13) for the lower 
face gives 
R 
l+u 
R 
y 
l = # 
The boundary conditions (A-12c,d) are difficult to 
satisfy. Using Hoff’s boundary conditions for uz,‘(x, ¥) 
and v,’(x, y) as in reference 6, and replacing (A-12c,d) 
with (A-13’), gives 
v = [(R/1 + ply atx = +al 
u’ = [R/(1 + pw) Ix y = +bf 


u'(x, y) = x + uz'(x, y) 
(A-13’) 


v(x, y) = + v,'(x, ¥) 


(A-14) 
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Nonsimilar Solutions of Equations (Continued from page 456) 





By. substituting Eq. (A-2) into Eqs. (11) and (12) 
and into the boundary conditions (14), (15) and 
(16) one obtains a set of equations which may be 
integrated from given initial conditions by a procedure 
similar to the one illustrated for the ideal gas case. 

For values of h/RT>) > 200, p/p; can be evaluated at 
a given station for specified values of pressure and 
temperature from a Mollier diagram—e.g., reference 
12. The integral of p/p, appearing in Fy, can then be 
evaluated as a subroutine in the program which solves 
the equations-~i.e., the second program. ** 
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Satellite Orientation by Inertial Techniques’ 
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Summary 


In this paper are derived the three-dimensional equations of 
motion applicable to accelerometers carried by satellite vehicles 
Utilizing the properties of the gravitational gradient, an instru- 
mentation system and an orientation procedure are proposed by 
which satellite attitude control and orbital parameter computa 
tion can be accomplished. An error analysis is included, and 
sensitivity and accuracy requirements are discussed 


Symbols 

Ay = acceleration sensed by accelerometer (vector ) 

A ipg = acceleration sensed by accelerometer located at 
Q with sensitive axis in the P direction 

G = gravitational acceleration vector 

AG = incremental gravitational acceleration vector 

g = amplitude of gravitation acceleration vector 

g = amplitude of gravitation acceleration vector at 
earth’s surface 

h = altitude above earth 

K = universal gravitation constant 

Me = earth’s mass 

mp = mass of particle P 

Rpg = distance vector from point Q to point P 

[R,.(a)| = nth rotation matrix through angle a 

R = distance from earth’s center to vehicle’s ¢.m 

rp = displacement from vehicle’s c.m. in P-direction 

ro = mean radius of the earth 

Y,, 1, Z; = orthogonal axis system, inertially fixed rota- 
tionally, and earth centered 

Xo, Yo, Z) = earth-centered orthogonal axis system with Z, 
through the vehicle’s c.m. 

Y,, ¥,,Z. = orthogonal axis system centered at vehicle’s 
c.m. 

13,2 = orthogonal unit vectors aligned with the Y, VY, Z 
axes respectively of the S-system 

ing = unit vector in direction from q to p 

w = scalar value of angular velocity of S-system with 
respect to U-system 

Ws; = component of ws in P-direction 

w = scalar value of angular velocity of 0-system with 
respect to l-system 

Wo} = component of w, in P-direction 

a, B, ¥ = angles of displacement of Xs, Ys, Zs; system from 
Xo, Yo, Zo system 

p = scalar value of distance of accelerometer from 
earth’s center 

Qpe = angular velocity vector of rotation of P-space 
from Q-space 

‘ = instrument error 

Je = error in measurement of quantity q 


Received April, 1960. Revised and received August 23, 1960. 
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Satellite Orientation by Inertial Techniques 


a DEVELOPMENTS and prospective applications 
of satellite vehicles in both scientific and military 
roles indicate the need for further study of the attitude- 
sensing problem under orbital conditions. Of the many 
techniques available utilizing radar, infrared, optical 
and inertial equipment, the utilization of the last per- 
haps is the most radically unconventional. It is the 
purpose of this paper to propose an inertial system 
based on gravitational gradient sensing and to discuss 
an orientation procedure and the accuracies required of 
the components as a function of the overall system ac- 
curacy desired. 

Conventional inertial systems using accelerometers 
experience difficulty in sensing the acceleration of a 
vehicle which approaches orbital conditions, due to the 
well-known inability of accelerometers to sense gravi- 
tational accelerations. Several authors! * ® have dis- 
cussed the possibility of sensing the gravitational 
gradient thru the use of accelerometers or stress-meas- 
uring devices. Since gravitational gradient sensing 
appears to be feasible, it will be discussed further on a 
three-dimensional basis including all vehicle motions. 
The basic sensing instrument is an accelerometer. A 
sequence of operations required to perform vehicle 
alignment is derived, and a computational method for 
determining vehicle altitude, orbital angular velocity, 
and angular acceleration is shown. 

For purposes of discussion, consider a vehicle in orbit 
about the earth which acts as a point gravitational 
source, and neglect the gravitational effects of other 
planets and the stars. The center of mass of the vehicle 
is the only point in the vehicle whose motion describes 
the orbit. Consequently, the total gravitational force 
on the vehicle results in equivalent inertial forces due to 
accelerations of this point. Since the gravitational 
force varies inversely as the square of the radial distance 
from the gravitational source, and since the centripetal 
force required for orbital motion varies directly as that 
radius, it is apparent that there must be forces exerted 
on other points of the vehicle structure in order to con- 
strain them to remain at a fixed displacement from the 
This may be illustrated by the simple 
The balance of forces at the 


center of mass. 
case of a circular orbit. 
¢.m. requires that 


—(Km,.m/r?) + mrw? = 0 


At a point in the vehicle structure displaced a distance 
Ar from the c.m., the force on the point is 


AF = [+(2Km,.m/r*) + mw*|Ar = [8Am.m/r*|Ar 








494 JOURNAL OF THE AEROSPACE SCIENCES 


Zo is 
Ru 
~~ 











AX, ge 
ig / An 

Xo SET(2) ‘Ss “ea 

Js 
Rig = RADIUS VECTOR—__ ted [Reas#)] Ray’ [Rev] Xo 
Ys Yo 
2s Zo 

Z, AZo 
Yi 
EARTH CENTER, 1,0 Yo 
= Xo. X 


v4 


Fic. 1. Reference axes for gravitation sensing 


which may be defined as the gravitational gradient 
force. This phenomenon may be used to advantage by 
sensing system from which the orbital conditions may be 
derived. 

To develop the three-dimensional equations of mo- 
tion, which include vehicle rotational motion, the 
Coriolis law is used to obtain accelerations at points 
displaced from the center of mass. Fig. 1 shows the 
geometry used. The coordinate systems are defined as 


follows: 

Xi, Yi, 2:1: Earth-centered inertial axes with X, 
parallel to the angular velocity vector de- 
scribing the orbital frequency; 

Xy, Yo, Zo: Earth-centered axes with Xo along X,, and 


Zy through the vehicle’s ¢c.m. 

Xs, Vs, Zs: Vehicle axes centered at vehicle’s c.m. 
displaced from Xo, Yo, Z) by angles a, 8, 
and y, as shown. 


For purposes of illustration, three sets of accelerome- 
ters are located each at a distance 7 from the c.m. along 
the X,, Y;, and Z, axes respectively. Each set, in turn, 
contains three accelerometers oriented with sensitive 
axes along Xs, Ys, and Zs. The following are defined: 


P : a point fixed in the S-space 
Qi: rotation of the 0-space from the 1-space 
Qos: rotation of the S-space from the 0-space. 


The analytical procedure used follows the methods of 
operation of the system. Consider that a vehicle has 
been placed into an orbit but that the characteristics of 
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the orbit were not known accurately ahead of time 
owing to uncertainties such as propulsion impulse and 
guidance inaccuracies, for example. Also, the vehicle 
may be rotating with respect to inertial space. Con 
sider that the task of the instrumentation system is to 
locate the local vertical, provide information by which 
the vehicle can orient its longitudinal axis perpendicu 
lar to the local vertical and in the orbital plane, and 
determine the characteristics of the orbit. The se 
quence used might be as follows: 

(1) Stop all vehicle rotation with respect to inertial 
space. 

(2) Align the vehicle so that its Xs axis lies in a 
plane perpendicular to the local vertical. 

(3) Align the vehicle so that its Zs axis lies along the 
local vertical. 

(4) While maintaining the last two conditions, align 
the Ys axis so that it lies in the orbital plane. 

(5) Compute the orbital characteristics. 

Reference to Fig. 1 shows that these conditions 
may also be represented respectively by: 


(1) Qs = —Qy 

(2) a =0 

(3) B =0 

(4) y =0 

(5) Computation of R, a, and aw» 


Generalized Vector Equations of Motion 
First develop the equations of motion with no restric 
tions placed on the vehicle rotation. Referring to Fig. | 


one may write: 
[Rip |i = [Ris 


Considering that the external forces acting at a point 
are composed of field and nonfield forces, one observes 
that 


1+ Reel, (1) 





[Rip]: = As + G + AG (2) 

[Rush = G (3) 

where G is the gravitational acceleration at the c.m., 

G + AG is the gravitational acceleration at P, and 

A, is the acceleration sensed by an accelerometer at P? 
due to forces exerted by the vehicle. 

[Rsp]: may be further developed by the Coriolis law as 


[Rsphi _ [Rsp lo + 2Qi0 X [Rsp lo a 
[Quo | x Rsp + Qw XK (Qi X Rsp) (4) 


Similarly, 


enh = [Rspls + 2Q0s X erie - 
[Qos]s X Rsp + Qos X (Qos X Rsp) (5) 


where [Rspo = snl + Qos x Resp (6) 


Noting that by definiton aks = [Rspls = 0, we 
may combine Eqs. (1) through (6) to obtain 


A, = (2Q1 + Qos) X (Qos X Rsp) + 
Bw X (Bw X Resp) + 
{ [2x0 ]o + [Qos]s} X Rsp — AG (7) 
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The three scalar equations for each accelerometer set may now be obtained by resolving Eq. (7) into the compo- 


nent parts. 


Matrix Transformation from 0-Space to S-space 
First it will be necessary to be able to resolve quantities defined in the (-space into the S-space. The coordinate 
transformation may be obtained by rotating 0-space coordinates through y, a, and 6 successively. Referring to Fig. 


1, one obtains: 





Ka l 0 0 ca 0 salfecy —-w OT X 

Ys| = {0 cB —s8 0 l O |} sy cy 0 | 

Zs 0 s8 c3 jL_—sa 0 ca || 0 0 1 |LZo J 
ca cy —casy sa Xo 
cB sy + sBsacy cB cy — sBsasy —ca sB |] Yo (8) 
sB sy — cBsacy sB cy + cB sa sy cB8ca |] Zo 


where the abbreviations s and c signify sine and cosine functions respectively. 
Evaluation of AG 
Elements of the matrix of Eq. (8) may be used to determine the components of AG in the S-space. Let 
[Rip| = p, IRis| = R, |Rsp| = 17 
then G = —(Km,/R*)(Ris/R) = —(Km,/R*)Ris and G+ AG = —(Km,/p*)(Ris + Resp) 
For accelerometer set (1), 
Ris + Rsp = (Rsa)i — (Reas@)j + (ReBca)k + ri 
p? = R°+ r+ 2rRsa 
I l 1/R' 
p : [R? + r? + 2rR sa]? . 


{1 + (*/R)? + (2r/R) sa}? 
| l | He ) 3r | 
- to Ee Dw oe 
p? ~—s«R X\R R 


Since R >>> 1, terms of orders of magnitude lower than r/R* may be neglected. Then one obtains 


1/p* = (1/R*)(1 — (8r/R) sa] 


By a series expansion, 


Hence 
= mM, ' R: R sa] sa}l R? R sa }(—casp) | Jj 
l 3r r 3r? 1 
R _ R sa }JcB ca |k + R - R sa if 


Noting that Gc = +Km,}sai — ca sBj + cB cak} 
and again neglecting the higher order term, one may write 
AG = Km,} (3r/R®)(sa)*i — (r/R*)i + (3r/R*)(sa)(—ca s8)j + (3r/R*)(sa ca cB)k} 

The value of AG at accelerometer set (2) is obtained similarly: 

AG = Km, (3r/R*)(—ca s8)sai + [(3r/R*)(—ca s8)? — (r/R*)]j + (3r/R*)(—ca s8)(cB ca)k} 
Again, in similar fashion for accelerometer set (3), 

AG = Km,} (3r/R*)(c8 ca)sai + (3r/R*)(cB ca)(—ca s8)j + [(37r/R*)(cB ca)? — (r/R?) |k} 
AG is now defined in the S-space for each accelerometer set. 


Resolution of Vector Equation Into Scalar Elements 


The remaining elements of Eq. (7) may be resolved as follows into elements defined in the S-space: 


cacy We, 
Qi = wo CBsy + sBsacwW |; Qs = | ws, 
sB sB — cB Sacy Ws, 
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Ca Cw Ws, 
[Q10 lo = wy cB sy + ssa cy ° [Qos ls = Ws, 
s8sy — cOsacy Ws. 
lr 
For clarity, Rsp is defined as Rsp = | 7% 
ts 
For set (1) rr =, r, = 7, = 0 
For set (2) n=! rr, =r, = 0 
For set (3) i a f, =e, = 


The cross-product operations can then be performed according to Eq. (7) for each accelerometer set. 


Complete Scalar Equations 

The total scalar equations may be written by combining the various elements. Since each accelerometer set is 
displaced a distance ‘“‘r’’ from the ¢.m., 7; = ’y =’: = /. 
The equations for the three sets follow: 


Set (1) 
Ag, = —2ww,,r(cB sp + sBsacy) — w,,?r — 2ww,7(s8 sp — CB sacy) — w,,°r — 
wr (s?B sta cty + s*y) — (38rKkm,/R*)s?a + (rAkm,/R*) 
A gy = —22apw,r cacy + w,,W5,% + wi*r(ca cy)(cB sy + sh sacy) + 
wr(sB s¥ — cBsacy) + a7 + (8rAm,/R*)(sa ca sB) 
Ay, = 2,7 CACY + W.,0,7% + wr Ca cy(sB sy — cB sacy) — 


ayr(cB sy + sB sa cy) — ar — (8rAm,/R*) sa ce ca 
Set (2) 


Agr, = 2wyw,,r(cB sp + SB sa cy) + ws, @,7 + wor (CB SP + SB sa cy)(ca cy) — 
ran(sB sy — cBsacy) — ra,, + (38rAm,/R*)ca s8 sa 


Ag, = 2wow.,r cacy — w,,°r — Zaw;,1(sB sp — cBsacy) — w,,"r — (9) 
wy*r[c?a c?y + (sBsy — cBsacy)*|] — (8r7Am,/R*)c?a s*B + (rKm,/R*) 
Ay,, = 2wow,r(cB sy + sBsacy) + w.,0..1 + wr(cB sy + sa sBcyp)(sB sy — cBsacy) + 
ray Cacy + @,r1 + (38rKm,/R*)c?a sB cB 
Set (3) 
Ay, = 20w,,7(SB sp — cCBsacy) + ws,0s,/ + wi r(cacy)(sB sy — cBsacy) + 
ran(cB sy + sBsacy) + a 1 — (8rKm,/R*) cB casa 
Ax, = 2wow.,r(sB sy — cB sacp) + w,.,Ws,1 + w*r(sB sp — cB sacy)(cB sp + sBsacy) — 
ron CACY — rwH,., + (38rKm,/R*) cB sB cea 
Ay, = 2,7 cacy — w,,°r — 2wyws 1 (C8 sy + sBsacy) — w, So 2 — 
wyr[e?a cty + (cBsW + sBsacy)*?] — (8rKm,/R*) 08 cea + (rKm,/R*) 
Since the amplitude and directional properties of the ac - io 7 a - 
gravitational field are implied in Eq. 9, this measure- 
ment system can provide the required intelligence for . ; , : . 
y I 8 For an inertially nonrotating vehicle, 2), = —Q);, so that 


vehicle orientation and orbital parameter computation, 


7 : : the right side of Eqs. (9) now contain only the gravita- 
which will be discussed next. 5 : ss 


tional terms. 
Vehicle Orientation Procedure The second step of the orientation procedure is to 
bring the Xs axis into a horizontal plane—i.e., one per- 
pendicular to the local vertical (Z) axis). In effect, this 
is done by reducing the angle a to zero. From Eqs. (9) 


The first step is to stop all rotational motion of the 
vehicle with respect to inertial space. This may be ac- 
complished with the use of rate-measuring instruments 


ve ; : one obtains 
such as rate gyros. With this done, one may write 


IR h = 0 Ay, = —(3rKm,/R*®) sa ca cB 
er Aag,, = —(8rKm,/R*) sa ca cB 
Thus, Eq. (1) becomes Ag, + Aaz, = —(38rKm,./R*) s 2a cB (10) 
[Rip]: = [Ris]: This value may be used to control motion about the J's 
be as ? axis. WhendA,., +A = 0, then o,. = o,, = 0, and 
and Eq. (7) simplifies to an te Rabie. SP, Sa 
1. (7) I a = 0° or 90°. The inertial stabilization is maintained 
A, = —AG about the Xs, and Yz axes, so that 
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= —w Ca cy, 
—wy(s8B sy — cB sa cy) (11) 


If a= O°. 
Ws, = —wty, ws, = —on SB sp 
If a = 90°, 
Ws, = OU, W:, = —w(sB sy — cB cp) 
The ambiguity in the value of a may be resolved by 


monitoring accelerometers .1 4, and «1 4,. 


If a = 90°, 


Ay, = wi'r(sB sp — cB ey)(cB sp + sB cy) — ras, 
= —ayw,,r(cB sy + sB cy) — ra, 
cl te = 2wyw (cp sy + s@ cy) + 


wy’r(cB sy + sB cy)(sB sy — cB cy) + asz,7 
= wow,r(cB sy + SB ch) + rex, 


Hence 4, + a, = 0. Ifa 0”, 
4, = wo’r(SB sp)(cB sy) — 

ran CY — rw, + (3rKm,/R*) sB cB 
Ag, = 2wyw,7 cB sw + 


wr CB sy SB SY + ray CW + ws,7 + 
(3rAm,/R*) sB cB 

= —wy"r(cB sW)(sB sW) + ray cy + @s,7 + 
(3rKms/R*) sB cB 


Ag, + Aa, = (87K m,/R*) 826 (12) 


vy 


Hence, if the total rotation about the X, axis with re- 
spect to inertial space is zero while the operation to cor- 
rect ais being performed, a monitor of 44, + A4,, will 
indicate whether a 0° or 90°. This logic operation 
can be used to obtain a = 0°. There is one further 
ambiguity, that a may be 0° or 180°, which apparently 
cannot be resolved by the inertial system. This will be 
discussed later. 

The third step of the orientation procedure is to re- 
duce 8 to zero, which will then align the Zs; and Z) axes. 
Eq. (12) provides the necessary information, but once 
again there is an ambiguity between 6 = 0° and 8 = 90° 
as well as integer multiples. The ambiguity between 
8 = 0 + amand B = (2/2) + an can be resolved by 
evaluating 44, + 44,. Noting that now w,, = @., 
5 = Ws, = a= 0,8 =OF 1, Tr 2 + mn, one may 


=“ 
obtain from Eq. (9) 
Ag, = —ray sy A 4y, = rao sy 


A Ay, + A 1, = 0) 

if 8 = 90°, and 

A ce = wo" sy cy 

Aay = 7 sy cy 

Aa, + Aay = wr s 2p (13) 
if B = 0°. 
Again the same logic can be used. If rotation about the 
Zs axis is kept inertially zero while the 8 correction is 
performed, a monitor of 44, + A ay, will indicate 


whether the value of 8 is 0° or 90°, so that simultaneous 
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monitoring of Eqs. (12) and (13) can be used to obtain 
B= 0+ nr. 

The fourth orientation step, alignment of the Vs and 
Vy axes, can now be performed using Eq. (13). Once 
again, the ambiguity in ¥ can be resolved by using 
other accelerometer outputs. Noting that now 


Ws = Ws = 0 a=£B=0 


Y= O0+n7, 7/2 +4 
and using -1,4,, as a monitor, if Y = 90° then 
A, = —@or 
and if y = O° then 
Ay, = 0 
The logic operation is again available to obtain y = 0. 


The fifth operation is to obtain information on the 


vehicle location. This information is now available 


from Eqs. (9) which have been reduced to: 


Aas, rKm,/R*® Ag, = 0 

Ag, = 0 An, —wy"r + (rKm,/R*) 

Ay, = 0 A 4,, = ra 

Aa © (14) 
Aay, = —Tro 

Ag, = —w’r — (2rKm,/R*) 


Assuming that the gravitational constant A and earth’s 
mass m, are known with adequate accuracy, the dis- 
tance R to the earth’s center may be computed from 
A, = rKm,/R*. However, the accuracy may be im- 


proved by using 


Ady, ~ Aas 3rKm,/R (15) 
wy) is available from 
14), = Peo 
or Aa, — Aa, = 2re (16) 
wy is available from 
Asay, + Aaz, = —2e0’r — (rKm,./R*) (17) 


Summary of Orientation and Computation Procedure 


(1) Reduce vehicle inertial rotation to zero. 


(2) From 44, + Aa,, = —(3rKm,/R*) s 2a cB 
(10) 
obtain a = 0° or 90° 
(3) Use 
Ag, + Aa = 9 when a = 90° 
= 3rKm,/S28 when a = 0° 


to resolve the ambiguity in a. Hence a = 0 + mr. 


(4) From Ay, + Aa,, = (3rKm,/R*) 526 (12) 
obtain B = 0° or 90°. 


(5) Use 
Aas, + Aay, 


0 when 6 = 90° 


0° 


= wor S2w when 8 = 


to resolve the ambiguity in 8. Hence 8 = 0 + mr. 
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1X 1025 components, angular errors ae, Be, and yY exist, and 
, r=|00FT. that these are the result of an error ¢ in each accelerome- 
ter, the following error equations may be derived from 
a 1x 1044 r=20FT. Egs. (9): 
i 7 _r=!0 FT. Ag, + Aa, = 0 = 2ay*r(Bepe — a) 
2 J—SFT. 24 HOUR (3rKm,/R®)2ae + 1.4e 
af j r=1FT. i (We + Beae) + (3rKm,/R*)2Be + 1.4e 
- "x Aa,, + Aa, = (1) = 2ay’r(We + Beae) + 
ml 1x 10°6- (3rKm,/R*)2aBe + 1.4¢ 
2 ; Eqs. (18) are set equal to zero since the system will try 
S 1 to satisfy such a condition. The term 1.4¢ is the root 
a IX 1°" sum square of two accelerometer errors. 
; As a first approximation, the angular errors may be 
] considered sufficiently small that the cross products of 
hike oe the errors are of lower order and may be neglected. For 
0. 1.0 10 100 1000 circular orbits, w)? = Am,/R*. The angular errors are 


ALTITUDE (10 FT.) 


Fic. 2. Gravitational gradient vs. altitude. 


(6) From 4,, + Aay, wy"r s2y (13) 
obtain y = 0° or 90°. 
(7) Use Aa, = 0 when y = 0° 
= —ayr when y = 90° 
to resolve the ambiguity in ¥. Hence y = 0 + ar. 
(8) Compute R from 
Ax, — Ag, = 3rKm,/R’® (15) 
(9) Compute w) from 
Ag, — Aay, = 2re (16) 
(10) Compute » from 
Aay, + Aa, = —2w0'r — (rKm,/R*) (17) 


Sensitivity and Accuracy Requirements 

It has been shown that alignment of the orbital ve- 
hicle and computation of the vehicle altitude depend 
on the ability to measure the acceleration, 3Am,r/R°, 
produced by the gravitational gradient. Although Eqs. 
(9) show that each accelerometer also measures quan- 
tities due to vehicle rotation, it can be shown that these 
are smaller than the gravitational gradient. Hence, for 
purposes of discussion of the sensitivity requirements, 
the latter effect is the dominant consideration. The 
gradient varies with the vehicle altitude and the separa- 
tion of the sensing element from the vehicle center of 
mass, as shown by Fig. 2. 

Two significant measurements to be obtained from 
the accelerometer system are the direction of gravity 
and the altitude of the vehicle. The effect of inac- 
curacies involved in these measurements can now be 
shown. 

Once the vehicle has completed the alignment proce- 


dure, ws = (0. Considering that due to inaccuracies of 


then 


0.53 e\ | 
ac = = ( 
3Km,/R?* \r 
0.707 € 
Be = % é (19) 
3Km,/R? \r 


212 (: 
aie 3Km,/ R3 ‘) 


The altitude measurement is obtained from Eq. 
(15): 


R? = 3rKm./(Aay, — Aaz,) 
The error is determined in similar fashion from Eqs. 
(9): 
A gn —_ As, ae | wor | (Bee =_ Qe)” = (We + aeBe)” } a 
(37rAm,/R*)Be? + 1.414} + (3rKm,/R*) 


Utilizing Eqs. (19), again neglecting the cross products 
of errors, and considering that the error ¢ is of lower 
order of magnitude than 3rKm,/R’*, the resultant error 


> 


in R is determined by 1/3 of the braced term in the pre- 


Re ATI (: 
R 3Km,/R* \r 
from which the altitude error is 


Ie | 0.471 * AC) (20) 
h 3Km,./R* h j\r a 


The errors in the remaining parameters defining the 
satellite orbit are obtained similarly: 


ceding equation: 


— ren. 1.115 (‘) (21) 
Angular velocity: -- = a 2 
6 “ @ 3Km,/R?* \r 


Velocity for circular orbit: 


Ve l a (‘) 
2 = : (22) 
J 3Km,/R? \r 


Angular acceleration: @, = 0.707 e¢/r (23) 
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SATELLITE ORIENTATION 


The errors described by Eqs. (19) through (23) are 
normalized with respect to the accelerometer error and 
the displacement of the accelerometer from the vehicle 
center of mass. They are shown as functions of alti- 
tude in Fig. 3. 

To examine typical errors for a 24-hour satellite ap- 
plication, consider a manned satellite where a 5-ft. 
separation is possible. Reference 11 indicates that 
accelerometers accurate to 10~!* g are feasible. Such 
conditions result in the following errors: 


a = 0.213 mr. w:/w = 0.04408% 
Be = (0.286 mr. Ve/V = 0.0632 
yp. = (),855 mir. 
h./h = 0.0223% @ = 1414 X 10-" 
rad./sec.” 
Conclusions 


Equations have been derived for attitude-sensing 
and orbital parameter computation in a satellite ve- 
hicle, and an orientation procedure has been proposed. 
A general accelerometer arrangement has been used 
for purposes of illustration. Eqs. (9) show the measure- 
ments obtained by orthogonal sets and can be used to 
determine an optimum arrangement if a specific meas- 
urement is desired. 

There is one drawback to the system proposed: the 
ambiguity between the angles of 0 and 180° in a, 8, or ¥ 
e.g., the system is unable to dis- 
However, a crude infrared 


cannot be resolved 
tinguish up from down. 
sensing device could be used to provide the required 
logic function for a and 8. Resolution of the am- 
biguity in y would require a time-history monitor. 

The analysis has shown the results that can be ob- 
tained through the use of accelerometers to measure 
the gravitational gradient. With accelerometers having 
a sensitivity range of 10~* to 10~'! ft/sec’, such as are 
feasible,'' '* and a 5-ft separation of the sensing ele- 
ment from the vehicle center of mass, the system can 
operate with good accuracy for satellites as far away 
as a 24-hour-orbit. 

In reference 2 a method is discussed for obtaining a 
desired accelerometer sensitivity over a wide range of 
values of gravitational gradient. Several sets, each 
having a different range, could be used to cover the ex- 
cursions of the gravitational gradient to be encountered 
by a vehicle due to the eccentricity of its orbit. In this 
manner, the system accuracy could be preserved for 
highly eccentric orbits. 

The first step of-the orientation procedure proposed, 
which is to reduce inertial rotation to zero, requires the 
use of high-quality rate-sensing instruments. The 
accuracy requirement is obtained from the condition 
that residual rotation due to rate errors must produce 
accelerations of lower magnitude than the gradient 
acceleration at the required orientation tolerance. By 
an error analysis similar to that discussed above, it can 
be shown, for example, that at 500 nm altitude, a rate 
error of 2 deg/hour will result in a, = 0.0047 rad. 
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This can be achieved by an instrument such as the MIG 
gyro. For comparable accuracies at higher altitudes, 
or for improved accuracies, more sophisticated rate 
sensors are required. For example, in a 24-hour satel- 
lite application, the same angular error requires a rate 
sensitivity of 0.014 deg/hour. A _ general class of 
solution would employ a rate-integrating gyro about 
which a closed-loop servo could realize measurements 
of the inherently more accurate rate information avail- 
able. 

The error analysis contained herein has been per- 
formed to establish sensing instrument requirements for 
gross system feasibility. Once these conditions are 
satisfied, the accuracies required of computer com- 
ponents used to process the instrument output data 
should be comparable to those of conventional systems. 
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On the Design of a Two-Dimensional 
Contracting Channel 


Y. S. Nanjunda Swamy * 

Research Assistant, Civil and Hydraulic Engineering Section 
Indian Institute of Science, Bangalore, India 

January 30, 1961 


1946, B. Thwaites! developed a method for designing an 


N 
I axisymmetric contracting channel with monotonic velocity 
distribution on its surface. In the present note, this method 
is applied to two-dimensional flows and a contracting channel 
is designed to have monotonic velocity distribution on its surface. 


The velocity potential #(x,y) in irrotational motion of an 


incompressible, inviscid fluid in two-dimensional flow satisfies 
the differential equation 
(0?&/Ox?) + (0?b/dy?) = 0 (1) 


axes ox, oy of the Cartesian system are 


axis of the contracting 


where the coordinate 
taken along, 
channel. Here, 
and x = a that of low velocity in the contracting channel. 


and perpendicular to the 
x = 0 will be taken as the end of high velocity 
The velocities “, v, in the directions of ox, oy are respectively 


given by 


u = db/dx, v = d/dy (2) 
Now, let us consider the flow taking place between two equi- 
potential planes situated at x = 0 and x 7, perpendicular 
to ox, with a potential difference of aor between them. Then, 
owing to symmetry of flow with respect to y about ox, the upper 


half plane y > O only need be considered. In this case, the 


solution 

* The 
Professor K. 
throughout the progress of this work. 


author is thankful to Professor N. 5. Govinda Rao and Assistant 


Seetharamaiah for their suggestions and encouragement 
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N 
ap. 
P(x,y) = ax + z ? sin px cosh py (3) 
p=1P 
satisfies Eq. (1) and the conditions laid above. Using Eq. (2) 
N 
uU(x,y) = a +r +» ap Cos px cosh py (4) 
p=1 
N 
v(x,y) = >. dp Sin px sinh py (5) 
p=1 


The stream function ¥(x,y) given by 


y 
V(x) = f udy (6) 
: ’ : 


i cos px sinh py ( 
> 


becomes 


Y(x,y) 


where p isa positive integer, ay’s are coefficients to be chosen 
to obtain monotonic velocity distribution on the surface of the 
contracting channel and N can take any positive integral value 


From Eq. (4) we have, 


y 
(*) --Fresnmompe 
OX J y={M p=1 


where M is a given number. If the coefficients a, be chosen 


such that 
cos x)*~! sin x, A>O (9) 


(0u/OX)y — M —A(1 4 


then, it can be shown for a given JN, that the first nonvanishing 


derivative of u, with respect to x at x w,y M is of an even 


order and is positive if a, > 0. Therefore u attains a minimum 
valueatx = 7,y = M. Attheendx 
derivative of u with respect to x at y 
and negative. Therefore, u has a maximum value at x = 0, 
y M. From Eq. (9) (0u/0x)y 
any other value of x in 0 < x < zw and, hence it follows that, 
u is monotonic decreasing in 0 < x < aw and y = 
M, < M. that v is also monotonic 


Therefore, the resultant velocity q f. q2)1/2 


= (), the first nonvanishing 
M is of an even order 


~ 0, and is negative for 


AM, where 
Then it can be shown 
(u? + vz)!’ is monotonic 
in the region considered. 

The following example worked out will illustrate the pro 
cedure for design. Let y = Y(x) represent the bounding stream- 
line which coincides with the surface of the contracting channel 
Let Y(O) = 1/3, V(r) = 1 = M. Let ap 0, p >4 so that 


N 3. Then from Eqs. (8) and (9), we get 

—} la, cosh 1 + (2a2 cosh 2)(2 cos x) + (8a3 cosh 3) X 

(4 cos? x — 1)} sinwxw = —A(1 + 2 cos x + cos? x) sinx (10) 

from which 

(a, cosh 1 — 3a3 cosh 3)/1 (4a2 cosh 2)/2 = (12a3 cosh 3)/1 
(11) 


From Eq. (11) on taking a; 1, we have, 


az = (2/5)(cosh 1/cosh 2) = 0.1641 

az; = (1/15)(cosh 1/cosh 3) = 0.01022 
a is determined by using the fact that the bounding streamline 
should pass through the points x = 0, y = 1/3 and x = -, 
y 1 and equating the values of the stream function at these 


two points. Thus, from Eq. (7) we have 


1/3a)9 + a; sinh (1/3) + (a2/2) sinh (2/3) + (a3/3) sinh 1 = 


a — a sinh 1 + (a2/2) sinh 2 — (a3/3) sinh 3 (12) 


With the values for a), ade, a3 given above ay works out to be 
1.9712. The value of the stream function along the bounding 
streamline will then be 1.0594. The ordinates of the bounding 
streamline are obtained as the solutions of the equation. 
3 
1.0594 = aoy + > > oP cos px sinh py (13) 
p=1P 
for allx in 0 < x < zw and with the values for do, ai, a2, a3 given 


above. 
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The method of solving Eq. (13) is to express it as a cubic in 
cos x and solve for cos x for every value of y on the bounding 
streamline. Hence x is known for each y on the bounding 
streamline. 

In the present example, the bounding streamline computed 
from Eq. (13) is shown plotted in Fig. 1 

The velocity distribution is given by Eqs. (4) 
streamline computed for the 


and ( 5) The 


velocities along the bounding 
example given above are shown plotted in Fig. 2 





Velocity 
distribution along the surface of contracting channel 


Fic. 1. (left): Solution of Eq. (13); Fic. 2 (right 


It must be mentioned here that the velocity distribution inside 
the contraction will be modified when parallel ends are attached 
and this effect may be determined by testing the contraction 


shape in an electric tank. 
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Fundamental Step Profiles in Thin-Airfoil 
Theory With Hydrofoil Application 


Andrew G. Fabula 
Naval Ordnance Test Station, Pasadena, Calif. 
September 15, 1960 


ip RECENT YEARS the linearized theory of steady, plane perfect 
fluid flow has been applied to various problems of cavitating 
and ventilated hydrofoils.' Such work 
essentially an extension of thin-airfoil theory by the complication 
In the course of some such work, a simple 


may be viewed as 


of free streamlines. 
fundamental-solution method of superposition of ‘‘step profiles’’ 
(of either camber or thickness) has been noticed which leads 
quite directly to lift, drag, and moment and which applies quite 
widely in thin-aero-hydrofoil theory. While final results should 


be identical to those obtained by usual techniques, the method's 


attractive simplicity recommends it for possible application t 
more complicated problems. 

For several basic classes of hydrofoil-with-cavity flows—e.g., 
the side-venting type? in Fig. 1, the physical plane (2 = x + cy), 
with foil and cavity linearized to the x axis, 0 < x < /, is mapped 
conformally onto the upper half of the “circle plane’ (¢ = — + in) 
with the wetted surface on ¢ = e’®, 0 < @ < wand the free stream- 
lines on ¢ = & |&) > 1 With U’., the freestream speed and 
U + iV the complex total velocity, then the conjugate-complex 
perturbation velocity and analytic function to be obtained is 
wz) =(U — U, — iV)/U The linearized Bernoulli equa 
ie, Pp — po = 
condition of constant 


tion without gravity —pU.*u, and the free- 


boundary cavity pressure 


streamline 
p. give 
p/2)U 7] 


t{ >1, uw(t,0) =u. = K/2, K =(po — ~) 


The linearized wetted-surface boundary condition on the upper 
(+) and lower (— ) wetted contours h,(x) is o(x,0+) = dh,/dx. 
Thus v is specified over the entire unit circle since analytic 
Aside 


-) = 


continuation into the lower half ¢ plane gives o(¢) = v(¢). 
from matters such as cavity termination, wake thickness and the 
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W (oo 14 
cu(x0+) df dx 
4 = 
23-| =( 2 
CAVITY WAKE I sa ae 
a 17 
s ; oa 
w (co)=C 
Fic. 1. (left); Fic. 2. (right). 


d(u-) relation, the problem of finding w(¢) — u, is similar to that 
for the camber line of thin-airfoil theory. Thus step-profile 
superposition will be illustrated by treating that well-known 
case. Let the camber-line contour be h(x), 0 < x < 1, so that 
the linearized, boundary-value problem on w(z) and w(¢) is as 
in Fig. 2. Let the step-profile camber line be h(x) = 0, 
O<x< », and A(x) = 6,x%,< «<1. Thus the step profile 
has v(e?) = 0 except at the only allowahle singular points of 6 
= m (the leading edge) and +0; = +ecvos~(2x, — 1). 
The function 


— 26 om” Mihai 2 
u(¢) = — -_ + + 
‘ rsin@ \¢ — e t— eh t+1 


is seen to give the step profile as follows. It has w(o) = 0 
and on the unit circle, except at the singularities, gives v = 0. 
The implied wetted-surface contour, by the Cauchy integral 


theorem, is 
h(x) = —Im a w(x )dz 


where C’ reaches from (0, 0) counterclockwise around the 
singularity slit to (x,0—) or to (x, 0+). Thus in the ¢ plane 
the only contributions to h,(*) occur as C’ passes e~*? and 
et which are seen to give, respectively, and as required 


— 2%6 dz\| 
—Im)| xi - eri (: = +§ 
7 sin 6 df]) ¢ = ¢* 


The leading-edge singularity, which makes the angle of attack 
zero, does not produce a step because dz/d¢ is zero for ¢ = e*”. 

The solution for arbitrary camber is to be obtained by re- 
placing 6 in w:(¢) above by (dh/dx,)dx;, and integrating along 
the chord. However this is not needed to obtain Cz, Cp, and 
Cy in the linearized theory where, after including the leading- 
edge singularity contribution to Cp, 


. L 

c= (p/2)U 7 = —2) $ wae = 4rb, 
; D f ‘ 

Cp = (p/2)U. = —Im ow dz = 0 

; M . f 

Cu = (p/2)U. = 2Rl Cc wedz = —4rbe 


where M is about the leading edge and positive clockwise; 
where C encircles the singularity slit; and where the expansion 
of w(z) for — © is 

wiz) = (a, + th)/z + (a2 + the)/s? + 


While Cp is trivially zero for fully wetted profiles, similar integrals 
are used for free-streamline flows with nonzero form drag. 
The step-profile camber-line solution gives (for the step at x) 


sites a2 + ih = +0 La (1 — 2x) 


d*- 2r 1— x 


; —16 


a Tr ib, 
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Thus the camber line A(x) gives 


1 an 
iow age uf ey 
0 dx V 1—-x 
tah | x 
Cu =-—2 2 (1 — 2x)dx 


0 dx V1 — <x 
and these are the same as obtained—e.g., from Fourier expansion 
of dh/dx* and other methods. 

Step-profile fundamental-solution superposition brings to 
such problems a little more of the directness seen in the thin- 
airfoil solution for a symmetric profile at zero angle of attack 
with thickness distribution ¢(x), 0 <x < 1, with 0) = #1) = 0, 


namely, 


1 ('da 1 
w(s) = - dx, 
2r 0 dx; = & 


This can be seen to be a particularly simple case of superposition 
of semi-infinite step profiles of thickness r = (dt/dx,)dx,, with 
the fundamental solution 
wz) = (7r/2nr)[1/(z — x)] 

The modification by conformal mapping to give free-streamline 
detachment from a cutoff trailing edge is straightforward and 
gives, after consideration of cavity termination, etc., the funda- 
mental step-profile solutions for hydrofoils with trailing-edge 


venting. 
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Experimental Data on Strain Accumulation 
Under Equivalent Thermal Cycling 


B. E. Gatewood,* A. P. Grothouse, and W. W. Von Hausen 
Air Force Institute of Technology, 

Wright-Patterson Air Force Base, Ohio 

September 16, 1960 


IT REFERENCE 1 calculated theoretical curves are shown for 
strain accumulation under thermal cycling for a two-element 
structure. Some recent experimental work?’ has been con- 
ducted on the three-element structure shown in Fig. 1. To 
avoid the effect of material property changes and the effect of 
temperature on strain gage and extensometers, the center upper 
element was heated instead of the test specimens. This pro- 
duced a thermal stress effect based on an equivalent temperature 
The strains were measured by extensometers on the three-test 
specimens while the stresses were obtained from strain gages 
on an elastic section of the upper elements (area ratio 3 to | 
in comparison to test specimens). 

The tests verified all the effects of thermal cycling predicted 
in reference 1: strain accumulation may occur under non- 
uniform temperature cycling, elastic shakedown may occur, 
strain divergence may occur, and allowable load-strain curves 
for the entire structure can be constructed. Fig. 2 shows the 
experimental results for an axial-load average stress of 43,800 
psi on the 2024-T3 clad aluminum-alloy test specimens under 


* Now, Professor, Dept. of Aeronautical Engineering, The Ohio State 
University, Columbus, Ohio. 
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formulation of their elastic behavior is imperative in structural where p and q are constants so chosen that ~ 
sails oe 2 se 4 . OL AT. mec eT] Ne 
design. The approximate solutions for several types of sections See — +1 + Ben = 0 and (1-2) © 1 9) ‘ 
have already been presented.!~€ Each of these has employed \ 
the Galerkin method. However, the Kantorovich method Phen Eq. (7) becomes 
employed in this paper is a generalization of the Galerkin d2u 5 q{5n? + 1) — 1 
method,’ which is very appropriate in dealing with the boundary- dz? = 20% 1 n)? + 122 “ur 
, , 2 2a*(1 — n)? 42? 
value problem peculiar to the geometrical configuration here 
considered. The solution to the problem will be approached 5 (¢\? g-a(p+2n-2)-2 (10 
‘ ; : ? 2 2n—2) —2 ) 
by use of the Prandtl stress function, applying the Kantorovich 2 \a N 
method. The use of the Prandtl stress function in pure torsion : 
problems has been outlined by Stanisié and Shirely.2 Therefore, Now let \ 
only the mathematical analysis peculiar to the Kantorovich c? = §/[2a%1 — m)?] and U1 + 1) = 1/4[q%(5n? + 1) — 1] , 
method and the stated problem will be presented in this paper. 
The Prandtl stress function ¥(x, y), is determined by finding (11) and 
the solution to Poisson’s differential equation Eq. (7) can then be written as how 
V(x, y) +2 =0 (1) , d%(uz—1/?) d(uz-"!2) ) ; ; Eq 
: 5 +3 — [c?2? + (1 + 1/2)?]us—? = bias 
under the condition dz? dz WA 
WL) = 0 (2) 2 (2) .-amten-m-12 (42) 

: , : P 2 \a whe 
where L is the boundary of the region ) under consideration _ ; : , The 
and V? is the Laplacian operator. For definiteness we shall con- Now letting / + 1/2 = » and employing the following two the 
sider that the region D is bounded by the lines x = Kj, and x = transformations successively : equi 
Ko, and by the curves y = g(x), y = A(x). The Prandtl stress y = uz? and 3? = (1/c?)Z? (13) 
function can be written in the form . Agi 

then Eq. (11) becomes 
n 
a a ee oy ‘ ey _, dy : 5 fa\2/Z \ —e t2n—-2) —1/2 : 
W(x, 9) = he xK«(%, Wf (x) (3) Z2-- 4 Z- — (Z24 yy = (‘ Sol\ 
K=1 dZ? dZ 2 r 
where the functions xx(x, y), (K = 1, 2, 3, ..., m) are equal (14) A : 
to zero on the contour Z with the possible exception of the line: a : ‘ ; 
‘ c poss € exc I ” be ; . The homogeneous part of Eq. (14) differs from Bessel’s equa- 
x = K, and x = Ke; and f(x) are functions satisfying Eq. (2) . . 2s me 
: ‘ ‘ : . tion only in the coefficient of y. However, it is well known that 
on the lines x = K,, and x = Ko. peice . ek : é , 
py : . e /“J,(iZ) is a real function of Z which is a solution to this E 
The Kantorovich form of the stated problem becomes s . 
homogeneous equation. It is denoted by the symbol det 
h(x) r 
f [V2y + 2)x«(x, y)dy = 0, (i= 1.2.2, ....:%) (4 x (Z/2)"+2™ I 
a(x) I(Z) = 5 
WZ) am (15) 
. ’ ' 7 mao @ir(y + m + 1) 
Let us now apply the Kantorovich method for the region, Fig. 1, 
bounded by the curves y = +ax”" and the lines x = K, and x = It follows that when » is not an integer, a fundamental solution wii 
Ke; where a, K,, and K2 are constants and nis any integer. For to the homogeneous part of Eq. (14) becomes 
the first approximation the solution is to be found in the form yn = AI(Z) + BI-,(Z) (16) 
of a polynominal of the second degree in y, and since it must 1 
vanish for y = -+tax", the Prandtl stress function is sought in where A and B are arbitrary constants. Now a _ particular Dia 
the form solution to Eq. (14) must be found. Letting Z = th, Eq. (14) Sep 
- es n is equivalent to : 
V(x, y) = (y*? — a*x*")fi(x) (5) 1 Hat 

, , , dy dy : P us F Rez 
corresponding to Eq. (3) with xx(x, vy) = (y? — a*x*") and n° th? II + (h? — v?)y = ah? (17) 
; i - 2 dh? dh ?’ 
fx(x) = fi(x). Eq. (5) becomes ner 

—_ where uw» = (3n + 1)/[2(1 — »)] and re 
[V2y(x, y) + 2](y? — a2x2")dy = 0 (6) - noe ee eee “’ 
—axn : ‘ 5 fqa\2/4\@to/2a—)) 
o= = (18) Wii 
Substituting Eq. (5) in Eq. (6) yields i al Ret 
d?f;(x) 5n df,\(x) 5n(2n — 1) a * 5x 2h It is easy to show that a particular integral of this equation, me 
dx? ” ‘ dx + 2x2 i: 292° 1 ii(x) = 2a? proceeding in ascending powers of 2 beginning with h#“+! is the 
es perl "a 
(7) Vp =a nee 
5 ‘ . (u + 1)? — »* Ela 
Employ the following transformation No 
hers 7 
- = gf -) = ug!? 8) : : . ae oe (19) ; 
x , fi(x) u (8) {(u + 1)? — vf f{(u + 3)? — v3} on 
(Rt 
- irae i es eo or employing previous substitutions, a particular solution /,(x) 
of Eq. (7) becomes 
J 5 " a, n— V/5n?+14+3 ; n+0J5n?+14+3 
; a Sr aaa , : 
8a(1 — n)? f 4(1 — n) 4(1 — n) 
iin~ = —™ [ et  _f. (20) | O 
m=0,, (—3n + 4m(1 — n) - VJ 5n? +147 » (30 + 4m = 2) + V/5n? +147 
4(1 — n) 4(1 — n) Ch 
Ae 
where I'(...) denotes the gamma function. Lai 
From Eq. (16) and prior substitutions, the homogeneous solution f;(x) becomes Se 
f(x) = Agi(x) + Byol 1] VP + BI V5/2 1-0) ga 1 , 
JAX) = Agi(x) J(x) = 44 -- - + lies , - gem p gli-sa)/2 (21) I 
V 5n?+1/[2(1—n)] \a(1 — n) — V 5n?+1/[2(1—%)] \a(1 — n) f 
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Geometry of the section. 


and a general solution, f\(x) is then given by the sum of the 
homogeneous solution Eq. (21) and the particular solution 
Eq. (20). Thus, Eq. (5) becomes 

(x, y) = (9? — ax" )\fi(x) = (y? — ax") [fi(x) + fo(x)] 

(y? — a*x*")[Agi(x) + Bgolx) + fp(x)] (22) 
where f,(x) and f;(x) are given by Egs. (20) and (21), respectively. 
The constants A and B in f;(x) must now be determined from 
the subsidiary boundary conditions, Eq. (2), and thus from the 
equations 
Aga) + Bg(a)+ f(a) = 0 and Ag,(b) + Bed) + f,(b) = 0 

(34) 
Solving these two simultaneous equations for A and B yields 


Sr(a)gi(b) — fp(b)gi(a) 


fy(b)gla) — fp(a)go(b) 


gi(b)g(a) — g:(a)go(b) 2i(b) goa) — gi(a)go(b) 


(35) 
Evidently the stress and displacement field are uniquely 
determined [see Eqs. (14-17)].° 
The torsional rigidity is given by 


b axn 
D = 2G f f Wx, y)dydx (36) 
a a—axn 


where G is the shear modulus. 
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On the Theory of Vortex Cancellation 


Clarence D. Cone, Jr. 
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September 20, 1960 


iy A RECENT PAPER! by A. Schaffer, the results of an experimental 
and theoretical study of the cancellation of a wing trailing 
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vortex by interception with another airfoil were published. 
The theoretical prediction of the angle of attack of the second 
airfoil necessary for maximum cancellation of the incident 
vortex was shown to agree fairly well with the experimental 
results. Examination of the theoretical development indicates, 
however, that a somewhat more exact treatment will produce 
even better agreement with the experimental data. 

The fundamental relation between the angles of attack of 
the leading (1) and trailing (2) airfoils is given in ref_rence 1 as 


ag = —a, + I'e:/7s' Vi (1) 


This supposes that a single resultant vortex of circulation —Is; 
will emanate from airfoil 2. However, the velocity distribution 
across the span of airfoil 2 due to the incident vortex produces a 
circulation distribution in which the circulation has a maximum 
value some distance from the median (center) section of the 
span. Therefore, the vortex sheet from the wing will consist 
of two bands of opposite vorticity which will form two discrete 
vortices upon coiling, one near the tip with I equal to the maxi- 
mum value on the wing, and the other further inboard with I 
of opposite direction and equal to the difference between the 


maximum- and median-section circulations on the wing. For 
maximum cancellation effect it is necessary that the incident 
and trailing vortices be superposed. Thus the outboard vortex 
of airfoil 2 is the important agent in cancellation of the incident 
vortex, primarily because the inboard vortex is relatively weak. 
Since the circulation of the outboard vortex must be equal to 
the maximum sectional circulation on the wing, the downwash 
velocity to be used in determining the proper a@ of the trailing 
airfoil must be the maximum, and this occurs at the inner 
periphery of the incident vortex core. Thus Eq. (1) should be 
modified to give 

ao = —a, + (T*1/27Vo)[1/a + 1(2s/ — a)] (2) 


where a is the radius of the incident vortex core of circulation I's;. 

The experimental determination of the actual radius of the 
incident vortex is complicated by the fact that the roll-up of the 
vortex sheet yields a vorticity distribution considerably different 
from a Rankine vortex.”* For an elliptical I distribution on 
the wing, in reference 3 it is shown that the sheet coils into 
This radius is con- 


semirotational cores of radius a@ = 2s/3. 
siderably larger than that for a Rankine vortex of equal kinetic 
energy, where a = 0.155s. Since the high velocities initially 
occurring near the center of the coiled-sheet core are rapidly 
reduced by viscosity, the maximum velocity will occur near 
the periphery or just inside it. Using the value a = 2s/3 in 
Eq. (2) with the model data of reference 1, the following is 


obtained for the three airfoil pairs: 


Pair No. A a) a (exp.) a (theory) 
1 4.00 10° —5.0° —5.9° 
2 1.94 14° —3.5° —3.6° 
3 1.16 18° —2.5° 3.2° 


The lack of agreement for the third pair is also observed in 
reference 1. It is possible that Eq. (5.2) of reference (1) is not 
sufficiently valid for predicting the relation between a; and Is; 
for such low-aspect-ratio rectangular wings; hence an over- 
estimation of I's; would overestimate the downwash angle of 
Eq. (2) and lead to overprediction of a,. By using the experi- 
mental values of a2 in Eq. (2), it is possible to estimate the 
radius within the core at which the maximum downwash velocity 
actually occurs. 

Finally, it should be noted that the linearized form of the 


downwash angle (I«:/2rVo)[1/a + 1/(2s’ — a)] must be re- 
placed by tan~!{ (I'*:/24Vo)[1/a + 1/(2s’ — a)]} for very large 


values of downwash. Then account should be taken of the 
inclination of the incident vortex with respect to the free-stream 
velocity both in the theory and in the experimental measure- 


ments. 
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An Experimental Study of Buckling of Thin- 
Walled, Pressurized, Conical Shells Under 
Compression and Compression-Bending 
Interaction 

Sq. Ldr. John K. Brown,* RCAF; Lt. Robert H. Rea,* USAF; and 
D. W. Breuer** 

Air Force Institute of Technology, Wright-Patterson AFB, Ohio 
October 5, 1960 


SYMBOLS 


E = Young's modulus 

M = bending moment 

P = recorded load at buckling minus load carried by internal pressure at 
critical section. 

y = radius of specimen 

S = section modulus at critical section 

t = wall thickness 

p = pressure parameter = (p/E)(r/t)? 

a@ = semivertex angle of cone 

oc = P/2rrt cosa 

om = M/S 

p = internal pressure 

«¢ = ocr/Et 

om = omr/Et 


onsen pressurized, conical structures appear in multi- 
stage rockets as well as in other pressure-vessel systems. 
There is very little published work, either theoretical or experi- 
mental, on the analysis of such structures. In order to gain a 
better understanding of their behavior, an experimental investiga- 
tion was conducted by the first two authors. It is hoped that the 
results of this study will be of benefit to other investigators. 

For the purpose of establishing a basis for comparison, one 
cylinder and three cones were tested. Fig. 1 shows the pertinent 
physical dimensions of all specimens. The wall material for all 
specimens was 17-7PH steel with E = 29.35 X 10° psi, ultimate 
strength about 200,000 psi and yield strength about 180,000 psi. 
The walls were mounted in 24ST-aluminum end plates with 
grooves cut for this purpose. The excess volume was filled with a 
low-melting-point alloy which served as the fastener between the 
end plates and the skin. The walls had two longitudinal seams 
diametrically opposed. These seams were one-half inch overlap, 
bonded and spot welded. 

All specimens were first tested in compression and then com- 
pression-bending interaction. The desired internal pressure was 
first developed and then external load was applied. For the 
compression-bending interaction tests the bending load was 
introduced after the internal pressure, and then the compression 
load was built up until buckling occurred. For the bending tests 
the wall seams were located at the neutral axis. 

The results of the compression tests are presented in Fig. 2. 
The ordinate on this figure is a stress parameter which is deter- 
mined by subtracting from the total external compressive load 
that part carried by the internal pressure. The radius and the 
area used in this parameter are taken at the cross section where 
buckling was initiated. The abscissa is an internal-pressure 
parameter. All symbols are defined in the list of symbols. It is 
to be noted that no data appear for specimen C; this is due to the 
fact that buckling of this specimen always resulted in permanent 
deformations, and a complete sequence of tests was not possible. 
For the other two cones and cylinder, the data shown were taken 
after three pretests. This was done because the first test always 
resulted in a buckling stress about 15 percent higher than 
succeeding tests. After this ‘“‘shake-down’’ phenomenon, repeat- 
able data could be taken. 


* Graduate students. 
** Associate Professor and Head of Department of Mechanics. 
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TABLE I. SPECIMEN PHYSICAL DIMENSICNS 


SPECIMEN re t 

CYLINDER |10" 5". 5" 0.005" 
CONE A io” +Ss° «4 0.005" 
CONE B oO s 3 0.005" 


CONE C Gc Ss 0.005" 


Fic. 1. Sketch of specimens. 


For all the data shown buckles were initiated at the minimum 
radius; however, if the pressure parameter is sufficiently small, 
buckling of the cones is not restricted to this location. Tests 
demonstrated this, but at these low pressures permanent creases 
remained after the removal of the load. It is to be noted that the 
data for the cylinder and cones exhibited a change in slope at 
about the same value of stress parameters but at different values 
of the pressure parameter. Although the slope is reduced at the 
higher pressure parameters, it is still well above the value of zero 
predicted by the theoretical curve. It can also be seen that the 
differences in cone geometry have apparently little effect on the 
compression data. The solid lines connecting the data points are 
simply best-fit experimental curves and do not represent theoreti- 
cal results. The data points from reference 6, shown for compari- 
son, are for a cone and cylinder with 707516 wall material and 
dimensions; Z = 20 in., 1 = 5in., r2 = 4 in., and ¢ = 0.004 in. 

The results of the interaction tests are presented in Fig. 3. 
The ordinate on this figure is the same as for Fig. 2, however, the 
abscissa is now a bending-stress parameter which does not include 
any pressure. The dashed portion of these curves show the 
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Fic. 2. Comparison of compression tests. 
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transition from pure compression to compression-bending inter- 
action. The spacing of the curves corresponds to the slopes of 
the compression curves on Fig. 2. The difference in slope of the 
curves as well as other differences in behavior between the curves 
for the two cones may be due to initial imperfections as well as end 
conditions, however, more testing is needed to clarify the situa- 


tion. 


REFERENCES 

1 Brown, J. K., and Rea, R. H., The Elastic Stability of Thin-Walled Pres 
suvized Conical Shells Under Compression and Compresston-Bending Inte? 
ictions, Masters Thesis, Air Force Institute of Technology, August 1960 

2 Donnell, L. H., and Wan, C., Effect of Imperfections on Buckling of Thin 
Cylinders and Columns Under Axial Compression, J. Appl. Mech., Vol. 17 
pp. 73-88, 1950 

3 Fung, Y. C., and Sechler, E, E., Buckling of Thin-Walled Circutar Cylinders 
Under Axial Compression and Internal Pressure, Journal of the Aeronautical 
Sciences, Vol. 24, No. 5, pp. 351-356, May 1957 

‘Harris, L. A., et al., The Stability of Thin-Walled Unstiffened Circular 
Cylinders Under Axial Compression Including the Effects of Internal Pressure 
Journal of the Aeronautical Sciences, Vol. 24. No. 8, pp. 587-596, August 
1957. 

» Lo, H., et al., Buckling of Thin-Walled Cylinder Under Axtal Compression 
and Internal Pressure, NACA, TN 2021, Washington, D. C., January 1950 

6 Lofblad, R. P., Elastic Stability of Thin-Walled Cylinders and Cones With 
Internal Pressure Under Axial Compression, TR No. 25-29 on Office of Naval 
Research Contract No. Nonr-1841(22), Massachusetts Institute of Tech 


nology, Cambridge, May 1959 


The Creep Deformation of Symmetrically 
Loaded Shells of Revolution 
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SYMBOLS 


‘1, % = principal radii of curvature 

€y, €@ = meridional and circumferential strains, respectively 
V, W = tangential and normal deformations, respectively 

t = time 

oy», 79 = meridional and circumferential stresses, respectively 


B,n = constants in creep law, Eq. (4) 
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é = effective creep rate, Eq. (4) 

a = effective stress, Eq. (4) 

P = pressure 

h = shell thickness 

a,b major and minor axes of e!lipse 
a a b 

Kk? (a? — 1) 


(1) ASSUMPTIONS 


it THE FOLLOWING ANALYSIS, the assumptions made are: 
(1.1) The sheil is thin and bending is neglected 
(1.2) Small-deflection theory is applicable 
(1.3) The shell is in second-stage creep so that a power-creep 
material law is suitable 
(1.4) The Soderberg combined-stress-creep relationships apply 
(1.5) The loading is time-independent 


(2) GENERAL ANALYSIS 
Using the equations for the strain of the middle surface (the 
notation of reference 1 is used), 

thé = OV/dOe — W (la) 
rep V cot g¢ — W (1b) 
Differentiate Eqs. (1) with respect to time, noting that the time 
derivatives of the radii of curvature are negligible as a conse 
quence of assumption (1.2). 
r,(O€,/Ot) (0?V/dtdg) — OW/dt (2a) 
ro(Oeg/Ot) = cot g(OV/dt) — OW/odt (2b) 

Eliminating (OW /ot) from Eqs. (2) gives 
(0? V/dtdg) — (OV /dt) cot ¢ 1(O€,/Ot) — ro O«eg/Ot) (3) 
Using the Soderberg relationships? for combined-stress creep 


together with a power-creep law 


‘js lo, — (1/2)o0] ¢/o é& [og — (1/2)e,]€/o (4) 
where 

o (og? + a9? — aaey) 

e Bo" 


B and n are material constants, and the dot indicates differentia 
tion with respect totime. Then Eq. (3) becomes 
(0?V/dtdg) — (OV /dt) cot ¢ 

Binloe — (1/2)o@] — re[og — (1/2)ey] jo"! (5) 
Now the right-hand-side of Eq. (5) is independent of time since 
the loading (and therefore the stresses) are assumed to be inde- 


pendent of time. Thus, Eq. (5) can be written 


0/d10V/O¢g — V cot ¢) Fig) 


where 
F(¢g) = Binloeg — (1/2)o9] — relog — (1/2)ey] 5 X 
lay? + a9? — ogo pias | 
Assume a solution for Eq. (6) of the form 
Vig, t) &(¢) T(t) (7) 
where 7(¢) must be a linear function of time, say of + « 
Without loss of generality, set c; = 1; further, since V(¢, 0) 0, 
it follows that c. = 0, so that 
T(t) =t (8) 


Substitution of the assumed solution leads to the following equa- 


tion for &( ¢) 
(db/dg) — cot g = F(¢) (9) 


for which the quadrature is 


&(¢) = sin of J [F(¢)/sin ¢elde + C} (10) 
Therefore the creep deformations are 
Vig, t) = tsin of J [Fl¢)/sin g]de + C} (11a) 
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W(¢, t) tcos gf J [F(¢)/sin g]dg + C} — 


Brot(og — (1/2)ey] (oy? + a9? — oyoe]-/2(11b) 
where Eq. (11b) has been obtained from Eqs. (11a), (2b), (4b). 
(3) CLAss OF SHELLS 
Consider now the class of .shells which satisfy the following 
conditions: (3.1) Closed shell of revolution. 
(3.2) Symmetrical about the equator. 
(3.3) Loaded with uniform pressure. 
From equilibrium considerations, the stresses are given by 
og = Pro/2h o9 = (Pre/2h)(2 — re/n) (12) 
Symmetry (3.2), requires that 
V(2x/2, t) = 0 (13) 
This condition may be used to evaluate the arbitrary constant 
C in Eqs. (11): 


= 
F(¢) ™/* F(¢) 

[Feat e- -f a= dg (14) 
sin ¢ e sing 


Therefore, the deformations are 


, on ee 7 
Vig,t) = —Btsin ¢ : dg (15a) 
¢ sing 
™/2 F(¢) 
W(¢, t) = —Btcos ¢ : dg + Btr(Pre/2h)” X 
yg sing 


[(r2/11) — 3/2] [(r2/r1)2 — 3(r2/n) + 3])-)/2  (15b) 


F(g) = 1e(Pre/2h)” [(re/n) — 1] X 
[(r2/r1)? — 3(re/r1) + 3] ~Y/2 


where 


(4) ILLUSTRATIVE EXAMPLE 


An example of the above class is the ellipsoidal shell with con- 
stant thickness for which the radii of curvature are 


m1 = da[K? sin? g + 1]~*/2 (16a) 
re = da[K? sin? g + 1]~”2 (16b) 
where a and } are major and minor radii, respectively, a = a/b 


and K? = (a? — 1). Then the deformations are 


w/2 
Vig, t) = —Bat(Pa/2h)"a"* K? sin ef g(¢) sin g de 
¢ 
(17a) 
W(¢, t) = —Bat(Pa/2h)"a"*[(1/2 — K? sin? ¢)g(¢) + 
‘9 


K? cos ef g(¢g) sin gdg (17b) 
¢ 


where 
e(g) = [1 — K* sin? g + K‘sin‘ g]~-/2/[1 + K? sin? g](+0/2 
Note that:for variable thickness, # must be included under the 


integral sign. 
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ANY FIRST-ORDER SLIP ANALYSES have been performed for 
boundary layers neglecting transverse curvature. It is 
the purpose of this note to examine the compressible boundary 
layer with transverse curvature in first-order slip flow. Correct 
to the first-order in slip, the continuum boundary-layer equa- 
tions with first-order slip and jump boundary conditions are 
applicable. The present analysis extends the continuum analysis 
of Probstein and Elliott! to include first-order slip. The Buse- 
mann and Crocco integrals are also discussed. No boundary- 
layer interaction effects are considered and only the zero pres- 
sure-gradient case is examined. f 
The axial symmetric laminar boundary-layer equations for a 
perfect gas with constant Prandtl number and specific heats in 
first-order slip flow are the same as the continuum equations 
. The con- 





given in reference 1 [see Eqs. (4)—(7) of reference 1 


tinuum boundary conditions, lim u(x, y) = ue(x)f, v(x, 0) 0, 
y—> a 
and lim 7(x, y) = T.(x), remain unchanged; the first-order 
yo 


slip and jump boundary conditions are?; 
u(x, 0) = al(ou/dy)y=0 + (3/4)[(u/pT)(OT/dy)]y=0 (1 
and 
T(x, 0) = Tu(x)t + bC(0T/dy)y=0 (2 
{being the mean free path near y = 0. 

The coefficients a and b are constants characterizing the type of 
reflections and amount of accommodation between the wall and 
the incident molecules. By an order of magnitude analysis, it 
may be shown that the magnitude of the second term, known as 
“thermal creep,” is of higher order than terms consistent with 
the boundary-layer approximation, particularly for highly cooled 
wall. In the following analysis, we shall therefore neglect this 


term. 


PARTICULAR INTEGRALS FOR Pr=1 
For Pr = 1 and isoenergetic free stream, the particular in- 
tegrals of Busemann (insulated wall) and Crocco (dp/dx = 0, 
and constant wall temperature) again may be derived, but are 
meaningful only if they satisfy boundary conditions correspond- 
ing to physically plausible situations. Consider the Busemann 
integral given by 
H =h + (1/2)u? = constant (3) 
At the outer edge, we must have H, = constant. At the wall, 
taking the wall derivative and multiplying by uw, Eq. (3) becomes 


PrxR(OT/Oy)y=0 + pwu(x, 0)(0u/dy)y=0 = 0 (4 


Maslen‘ pointed out that Eq. (4) with Pr = 1 is exactly the 
condition for zero heat transfer in slip flow when the work done 
by sliding friction is considered. Thus the Busemann integral is 
valid in first-order slip flow under the same conditions as in con- 
tinuum flow. The Crocco integral is given by 


H= Bu+cC (5 


Using the boundary conditions [Eqs. (1) and (2), and H, = const.|, 
we obtain from Eq. (5), respectively, 

hy + btB(Ou/dy)y=0 — Bat(dou/dy)y=0 + O(#) = C 

and Buet+tC=H, 


We see that in general the wall boundary conditions, Eqs. (1) 
and (2), cannot be satisfied since B and C are constants. It is 
observed, however, that for first-order slip flow and with a = } 
the constants B and C can be determined and are equal to those 
determined for the continuum case. The Crocco integral is thus 
valid only if the condition a = b is added to the usual continuum 
restrictions. 


THE CASE OF ZERO PRESSURE GRADIENT, UNIFORM STREAM, 
AND CONSTANT WALL TEMPERATURE 


Considering now the case of uniform stream and dp/dx = 0, the 
t+ Certain incompressible pressure-gradient cases can be analyzed simi- 

larly.? 
t Subscripts e and w correspond to edge and wall conditions, respectively. 
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modified Mangler (, 7) and the Howarth (%, Y) transformations 
may be applied to the equations of motion and the boundary 
conditions. Assuming a linear viscosity law and by introducing 


y VW (Coveted) Lf (%, Y), T/T. = »\(&%, Y), 9 = V (tte/Ceve) Y, 
and & VJ (Covek ue( L/ro2) cos a the Probstein-Elliott equations 
ire obtained [see Eqs. (31) and (32) of reference 1). 

The variable — is proportional to the ratio of boundary-layer 
thickness to body radius and is defined as the transverse-curva- 
The modified 


ture parameter. (For Mangler bodies, § —~ 0 


boundary conditions for constant wall temperature become 


Of Yo Pu U, O*/ 
£0 0, at 
On/ =0 i Vous On? /=0 
as 
1 ( 
On / n= 
ry p Uu, Or r. 
rn(E, O) = Ay + OE , ACE, Ll 
| Vo, ¥ \On/n=0 


For the same class of ‘‘admissible’’ bodies as in the continuum 
case, the solution for first-order slip may be written in terms of 
double asymptotic expansions in £ and ¢, where e, the slip pa- 


rameter is defined as al(m/L)V ue/CvX; ie., 


A(E, €, n) [E‘Xso(n) + Edu (n)] (9 


In the above, the powers of e greater than one are excluded 
since we are only considering first-order slip effects. Expanding 


the boundary conditions, Eqs. (6) and (7), in terms of Eqs. (8 


and (9), the boundary conditions become *** 
fio(O 0, fi() = 0 for i>o0 
0 0, f’i1(0) (1/Aw)f” io(O for i20 
(10 
s fsa 0 for i> i. | 
f-.( ) 8] for 170 
and 
Aw (O Aw, Aio(O O te «2 i, 
di (0) (b/a)(1/rAw Ai 6 (0 for i>0 
ei (11 
A = | Viol =@ for 21, 
Au(e) =O for 7 >0 


Substituting Eqs. (8) and (9) into the Probstein-Elliott equa- 
tions of motion and equating like powers of — and =< to zero, 
two infinite families of ordinary differential equations result 
The zero-order equations are the usual compressible flat plate 
equations. The first-order slip equation from the momentum 
equation is essentially the system obtained by Lin and Schaaf? 


for first-order slip on a flat plate. Its solution is simply fi 


, 


L/w) foo 
equation from the energy equation is 


, where foo is the Blausis function. The first-order slip 


(1/Pr)N"o1 + (1/2) foor’cr + (1/2)forr’oo 4+ 
Ay — IMA uf’o = 0 (12 


with boundary conditions given by Eq. (11 The solution for 


Eq. (12) with Pr = lis 
dor (1/rAw){A’o0 + [1 — (b/a)] [1 + M2 (y — 1)/2 — Awl X 
[f”o0(0)](f’0 — 1)$ (13) 


The first-order transverse-curvature equations in both the mo 
mentum and energy equations are the systems examined by 
Probstein and Elliott in reference 1 

The terms f/f}; and \, give the combined effects of first-order 
transverse curvature and first-order slip on the velocity and 
temperature profiles, respectively. From the momentum and 
energy equations, f;; and \,, are given by 


*** Primes denote differentiation 
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(1/2 — x)(f'oof’n — fuf”’o + f'af'1 — fof’) — 
(1/2) (fof +f’ ufoo + faf’10 + fiof’ — f'n 
” a) : 
21 foo Aoidn + f dn 14) 
J0 0 
and 
1/2 — x)(f’owrAn + for — (1/2)(foor’n + for 
1 — K)(fwrA‘an + fir’ — (1/Pr)y 
q n n 4 
(2 Pr) Xe Aor dn + A‘q Aov dn | + 2(4 ] x 
0 0 


n aa 
M,2 ff 0 + f'n" + fw? dn + 26" of” | Aoodn 
JV J 0 


(15 

The boundary conditions can be obtained from Eqs. (10) and 
(11). For ‘‘admissible’’ bodies, « is a constant.! 

In terms of the Probstein-Elliott continuum solution, /\{; may 


be written as 


Su (1/Aw) f'10 + Of 16 
where “f,; satisfies 
(1/2 — x)(f’o MP f'n — fu fo 1/2)(fuf’o 4 iz 
*n X "7 
“ AD gl Aor dn f 17 
J 0 d J 
with f,(0) = 0, DF’ (O 0 and DF, 0. The 
solution of Eq. (17) for Pr = 1 is 
Pee 
u = fo yA lidn + (1/Aw) [1 + 7.2 (y — 1)/2 r 
\ J0 
*” 
[1 — (b/a)]f 0 (Js J, )dn 
0 


” bd] FF ae A uy) a ) 
in | zr f’ wo? foolidn dndn (18 
0 0 f'00* ly? /7 0 f 


The functions /;;, J;, J2 and their integrals are obtained in refer 
ence 1 for the cone and cylinder with Pr 1. The constant A 
is evaluated using the boundary condition, f’ 0) 
For a cone, we have evaluated this solution.? 

The solution for \,; is not undertaken; however, we have al- 
b and Pr l 


Thus the solution for Eq. (15) 


ready shown that for the special case when a 
the Crocco integral is valid 
fora bis simply 
An [11+ (y¥ -—1)/2 M2 —d fy 
(y — 1) M2A(f'0 f'n + fats (19 

where “f;; is given by Eqs. (16) and (18) when a b 

The results of practical importance are the skin-friction 
coefficient and heat transfer to the wall; these are given by 

: Ou 
Cr = Mu [(1/2) pete” ] 

Oy/ y=0 
Cove 4 ” 


ueX L i=0 
and 


g = K(OT/dy)y=0 + pou (x, 0)(0u/dy)y=0 


Yo Pu 
Ue Ser. ) % [E*r’ (O) 4 E'eX i(0)] +4 
Lm Vous { i=0 


M Ue” trey cn” 
Eel f” o(O)f "(0 


Aw i=0 


+ f"-1,0 Of") +...4 


For the special case when a = band Pr = 1, it can easily be shown 
that the skin friction and heat transfer are directly related by 
the usual continuum Reynolds analogy. 

The skin-friction coefficient through the first-order in trans- 


verse curvature in continuum flow (« = 0) has been completely 
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determined by Probstein and Elliott! for cones (x = 2/3) and 
cylinders (x = 0) with Pr = 1. Using these results, our calcula- 


tion for a slender cone with Pr = 1 in first-order slip flow gives 
Cy/Cyy, = 1 + £[0.5166 + 0.9126d, + 0.1206(y — 1)Me?] — 
&€} 0.7996 + (0.3033/A~)[1 + (y — 1)/2M2 — Av] X 
[1 — (b/a)}} +... (22) 


For Pr = 1 anda = 5b, the Reynolds analogy gives for the heat 
transfer to the wall 
(Cy/Cy 1a =b = (g/gu)a=b = 1 + 
£[0.5166 + 0.9126 r» + 0.1206(y — 1)M.2] — 0.7996£e 

In the above, Cy yy and q,, are the continuum Mangler quantities. 

Without pressure gradient, the first-order effect of slip on skin 
friction enters first in the term of order £e, which is essentially 
the Knudsen number based on local body radius. For the case 
when a = 3, this is also true of heat transfer. Through the first- 
order in slip alone (£ = 0), these results agree with all first-order 
flat plate investigations, however, there is a slight contradiction 
to Maslen,4 who concluded that there was no first-order slip 
effects on heat transfer even with a # b. According to the 
present analysis, from Eqs. (13) and (21), the first-order slip effect 
on heat transfer without transverse curvature is 


e(1 — b/a)f"00(0)(9y,/Aw) 


which vanishes only if a = b. 
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HE SOLUTION of thin elastic shells of revolution based on 

Love’s first approximation! involves two differential equa- 
tions of second order generally known as Love-Meissner equa- 
tions. Since the inconsistencies of Love’s equations of equilib- 
rium as well as his approximate formulation of the stress-strain 
relations have been critically re-examined by numerous authors, 
an attempt is made to obtain a new set of governing differential 
eguations on the basis of a better second approximation. Such 
equations are presented here for thin elastic shells of revolution 
of uniform thickness subjected to axisymmetrical loads. Unless 
otherwise noted, the notations and sign conventions follow those 
used by Timoshenko and Woinowsky-kKrieger.*? Thus, when 
additional terms due to the second approximation are omitted, 
the equations will reduce to the Love-Meissner equations given 
in reference 3. 

For an element of a shell of revolution shown in Fig. 1, the 
consideration of extension and rotation of the middle surface 
leads to curvature terms in the equations of normal forces, 
strain terms in the equations for moments, and rotational terms 
in the equilibrium equations which are generally omitted in 
Love’s first approximation. Thus, the normal forces and 
moments can be expressed in terms of strains and changes in 
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Fic. 1. Equilibrium of an undeformed shell element. 


curvature as follows: 
No = [Eh/(1 — v?)], eg + veg + (h?/12r) X 
[(1/n)(dV/d@ — €g) + (v/r2)(V cot @ — €)|{ (la 
No = [Eh/(1 — v?)| beg + veg + (h?/127r,) X 
[(1/r2)(V cot @ — €) + (v/n)(dV/dd — eg)|} (1b 
My = [|E£h8/12(1 — v?)][(1/n)(dV/d@ — eg) + 
(v/r2)(V cot d — €9) + (1/re)(eg + veg) | (le 
Mog = [Eh§/12(1 — v?)|[(1/re)(V cot @ — eg) + 
(v/r)(dV/dp — eg) + (1/ri)(eg + veg] (1d) 
in which V is the angle through which the tangent to a meridian 
rotates during the deformation. The equations of equilibrium 
are obtained through the summation of forces and moments in 
the directions of the coordinate axes as follows: 


(d/do)(Ngrz sin @) — 1Ne(cos @ — V sin d) — 
Qgr2(1 + dV/dd) sing = 0 (2a 


(d/d)(Qgr2 sin @) + rNg(1 + dV/d¢) sin @ + 
rNg(sin @¢ + Vecos ¢) + Znr sing = 0 (2b 


(d/dod)(Mgrz sin ¢) — 1, Me(cos dé — V sin d) — 
Qenr2 sind = 0 (2c 
In order to simplify the solution of the problem, however, it is 
necessary to determine the relative importance of the terms 
introduced by the consideration of extension and rotation of the 
middle surface of the shell. It is noted that the curvature terms 
in Eqs. (la) and (1b) for the normal forces are multiplied by a 
factor of the order of 42/127? where h is the shell thickness and 
If the ratio of h/r ?s smaller 
than 1/10, a reasonable upper limit for the application of thin 
Thus, 


ry is the least radius of curvature. 


shell theory, this factor becomes smaller than 1/1,200. 
these terms in Eqs. (la) and (1b) are smaller than the remaining 
terms in the same equations. Conceivably, if V and its deriva 
tive are both large, the introduction of the curvature terms be 
comes significant, but this possibility appears to be remote. It 
is justified therefore to omit the curvature terms from Eqs. (la 
and (1b). On the other hand, the strain terms in Eqs. (lc) and 
(1d) for the moments have the same order of magnitude as other 
terms in these equations and hence should be retained. 

In Eqs. (2), the rotation V contributes terms involving prod- 
ucts of strain and derivatives of strain only." Since the dis 
placements are small, these terms are of a higher order and in 
general need not be retained. The consideration of such terms 
is necessary only in the investigation of large deflection of shallow 
shells with small ratio of rise to radius of curvature, and in the 


study of buckling of shells. 
After omitting all higher order terms consistently from Eqs. 
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1) and (2), a set of equations representing a second approxima- 

tion of Love’s theory is obtained. It is noted that the only 
difierence between the two approximations lies in the stress- 
strain relations for the moments in which the strain terms are 
included for the second approximation. 

On the basis of such simplifications, the stress-strain relations 
and the equilibrium equations may be reduced to two second- 
order differential equations in a manner similar to those based 
or Love’s first approximat on obtained by Meissner. Using 
Vand U = rnQ¢ as independent variables, and considering Z as 
the intensity of pressure normal to the middle surface, the fol- 


lowing results are obtained: 


L(U) + (vb/n)U = EhV + ® (3a 
L(V) — (v/n)V = U/D — H(U) + % (3b 


in which 


; ry dd? r 1 d [rz 4 ry a...) 
; 1 ot = 
r2 dd? r, Ldd \n ra we do 


r, cot d 
( 


Pet) 
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H 1 Yov (’ l ) d?(. ) ie 1 lov d 
rn? LEh \re ry dq? r, | Eh do x 
{1/1 1\) (2 — ot d fl l 
B¢-9)-g9@- 
(nn \re rJ \ Eh re r 
G22) (1 aces ( 1 
do r.Eh "2 ’ 
rm 1d (cot ¢ (: 1\) 
d —_ { 
Ehr,do@ | fr r2 nj) 
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fib [ Z rir COS @ sin @ dd 
J0 
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® ¢ (. se =) (Dp _ Zre? rs : ss Cc QP 
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* ( + ) — vZrs + Zr2 . x 
sin? @ \r ‘ Eh 
( 1 ) vry1 d ( 1 
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Consequently, the shear, normal forces and moments can be 


expressed in terms of variables V and U as follows: 


Qe = U/n (4a) 
No —(U/r2) cot @ — Zr2/2 (4b 
No —(1/r,;)(dU/dd) + Zro2/n — Zre (4c 
Me —D\|(1/n)(dV/dd) + (vV/rz) cect @ + 


(1/Eh)(1/rz — 1/m)(No — vNe)] (4d 


Me —DI\(V/rz) cot @ + (v/n)(dV/dod) — 

(1/Eh)(1/rz — 1/r)(Ne — vNg)] (4e) 
When the shape of a shell of revolution is defined by its meridian 
curve, Eq. (3) can be solved either analytically or numerically 
Hence, the shear, normal forces, and moments can be obtained 


by Eqs. (4 
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Note on Airfoil Theory in Hydromagnetics* 


G. S. S. Ludford 
Div. of Applied Mathematics, Brown University, Providence, R.I. 
October 10, 1960 


N THE SEARS-RESLER THEORY! of thin airfoils in fluids of in- 

finite electrical conductivity, the following potential problem 
arises when the applied magnetic field is perpendicular to the 
free stream. Determine a vortex distribution k(x) on the in- 
terval —b < x < b, whose potential ¢ satisfies the condition 


(Og/Ox)y — m(0¢/OY)y —mUY'(x (1 
and has finite derivatives at the trailing edge x = b} (Kutta- 
Joukowski condition). Here y = Y(x) is the camber line of the 


airfoil; U is the undisturbed flow velocity, and m (magnetic 
Mach number) its ratio to the Alfvén velocity determined by the 
applied magnetic field. 

In their paper, Sears and Resler use Glauert’s series 


> 


(O0¢/Ox)yet+o = R(x)/2 By tan 6/2 + > B,, sin n6 
n=1 
x = beos#A0 S<A<-w) (2) 
Then (O¢/Oy)y-0 = —Bo — p> B, cos né 
n=1 


and the condition in Eq. (1) becomes 


x 


> B, (sin n@ + m cos nd 


n=1 


—mUY'(x B,(tan 6/2 + m) + 
It is now a question of expanding Y’(b cos @) in the series appear- 
ing on the right side of this equation. Even if this can be done 
for (say) camber lines with continuous slope, there is no simple 
rule for determining the B,, since the functions sin »@ + m cos 
n6 are not orthogonal in the interval 0 <6 < 7. Nor is the in 
verse problem of finding the camber line, given the lift loading 


2 


I(x) = [(1 + m?)/m?] pU R(x (3) 


solved satisfactorily in this way. For although the B, can be 
determined explicitly from (2), it is not immediate'y apparent 
what functions k(x) are appropriate. Moreover, to determine 
the change in /(x) due to a given change in angle of attack of the 
airfoil so found, the direct problem must be solved for the case 
V’(x) = const. 

In fact, we shall see that the Glauert series is not appropriate 
for regular camber lines, although all that is needed is a simple 
modification of it. In general, the corresponding k(x) behave 
like (6 + x)~ at the leading edge and like (6 — x)? at the trailing 
edge : 

The starting point is the solution for yg, due to Rott and 
Cheng,’ which is given in the Sears-Resler paper’: 


o¢ 0¢ : mU (: _ sy 
—1 = 4 7 
Ox oy rV1 +m? \2 t b 
f (; 4 ‘ Y'(s J 
af b § sS—gs 


where z x + ty and 6 (1/r) tan-!' m. SinceO<m< @, 


we have0 <8< 1/2. Asy— +0 the integral tends to 


b (; +. ‘ V’'(s)ds ( + x\8 - 
+ m1 Y’(x 
a =" b — S, S$ = Zz 4 —_— £ 


since it is proportional to the complex velocity of a vortex dis- 


tribution with density [(6 + x)/(b — x)]® Y'(a The Cauchy 
principal value is implied. Also the factor [(z — b)/(s + b \8 
tends to [(b — x)/(b + x)]® e8*. Hence 
0¢ 2mU sin Br (b — x\? 
k(x) = 2 = 7 x 
Ox/y = +0 V1l+ mm? 7 b+x 
b b+s B V’(s)ds ad 
+ cos Br Y'(x (4) 
= b-—s s—x 


* This research has been supported by the Office of Naval Research 
under Contract No. 562(07) with Brown University 
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From this formula we can determine the character k(x) when 
With the possible exception of the lead’ng and 


Y’(x) is regular. 
To determine its behavior 


trailing edges it is also clearly regular. 
at these points, we need only consider the case of a flat plate 
Y'(x) = 1. The integral in Eq. (4) is then [s = b(2¢ — 1)] 


>} . t B dt Qh 
20 — x 
J0 1—t/ 2bt —- (6+ x) +x 
Br - 2h 
3 Ft i, 1 +262: 
sin Br b6+x 


where the mean of the values at x + O07 must be taken for the 


hypergeometric function F. Now 
F ~ (cos 78/8) [2b/(6 — x)|? — 1/8 for b — x small 
F ~ —(1/8)((6 + x)/2b| for 6 + x small 


Thus at the trailing edge, x = b, k(x) tends to zero, 


k(x) ~ —(2mU/V/1 + m2) [(b — x)/2b]8 


while at the leading edge, x = —d, it has the singularity 


R(x) ~ —(2mU V1 + m2) [2b/(b + x)]8 
This particular k(x) is of some importance; as mentioned above, 
it gives the change in lift loading due to change in angle of attack. 
We turn now to an effective representation for k(x) in general. 
Note that in the case just discussed 
[(6 — x)/(b + x)]! B I(x) 


+h as the lift loading itself in the 


has the same behavior at x = 
1/2), for which Eq. (2) is 


nonmagnetic theory (m = ©, B 
appropriate. Moreover, the terms tan 8/2 and sin n@ in Eq. (2) 
are separately produced by the integral term in Eq. (4)—s 
1/2—when 1 and cos n@ are substituted for Y’; this is the essence 
These two facts suggest using [(b — x 


of Glauert’s method. 
(b + x)]/2-8 Y"(x) here 

Thus the rule is to form the Fourier expansion 
— U(tan 6/2)!-28 Y’(b cos 0) = By + 


> B, cos nd (0 <9 < x) 
n=1 


then [see Eq. (4)], 


k(x) + (2mU/V/1 + m?) cosBr Y'(x) =(2m V1 + m?) X 


x 


sin Br(cot 6/2)'-28 [By tan 0/2 4 > B, sin n6} 


n=1 


and the lift loading is given by (3). 
It follows that the ideal angle of attack a is given by 


3 6\'? 
f tan [Y’(b cos 0) — ald@ = 0 
0 & 


Satellite Orientation by Inertial Techniques 


-JUNE 1961 


that is 


sin Br g\1—28 e 
tan 5 Y’(b cos 0)d@ 
us 0 2 


However, the lift Z and pitching moment J/ are not readily com- 
puted from Eq. (5) since the integrals arising from the series are 


not simple. Reverting to Eq. (4), we find from (3), 


l+m 
L= pl k(x)dx 
m? mi 
. } / 
V1 4A gy? _ | sin Ba iad b 
—2 pU? 
m —~§ 
b 
f y AS c : : 
ab (b 4 5 B-1 
oU2 | 2bB x 
m -_»p (b —s)! 


2b 
Fi 1,1 — 8, 2; 
bites 


Similarly 
1 -++ m? ; b 
pl xk(x)dx 
m? —b 


e/ 


b \ B-1 
2 (bh = .§ 
1 m? : t 
+m" oU?| bie x 
» (b— 8 
2b 
na Ss 
2h : 
23. > V'"(s)ds + 
b+sJ f 
b 
cos Br [ x Y'(x)dx 
e b 


Simplification of our results is contained in ref. 5. 


Y’(s)ds + cos Br} Y'(b _ 
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STANDARD-THOMSON CORPORATION 
STANLEY AVIATION CORPORATION 
THERM, INCORPORATED 

THIOKOL CHEMICAL CORPORATION 
THOMPSON RAMO WOOLDRIDGE INC. 
TRANS WORLD AIRLINES, INC. 


UNION CARBIDE METALS COMPANY, DIVISION OF } i 
UNION CARBIDE CORPORATION ef 


UNITED AIR LINES, INC. 

UNITED AIRCRAFT, CORPORATION 
UNITED STATES AVIATION UNDERWRITERS, INC. | 
VEHICLE RESEARCH CORPORATION 

WESTERN GEAR CORPORATION 


WESTINGHOUSE ELECTRIC CORPORATION 








